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Preface 


This report has been prepared in two volumes, each of which is a separate 
document. Volume 1 is in the form of the usual final report. It contains, a 
summary of the theoretical derivations, the required analytical boundary 
value solutions, and a numerical analysis of the solutions, as well as con- 
clusions and recommendations for further work. It includes all the equations 
needed to evaluate any of the boundary value solutions except those equations 
which apply strictly to two-body motion and can be found in most standard 
astrodynamics or celestial mechanics textbooks. 

The actual derivations of the second order asymptotic solutions are long and 
involved. These derivations have been compiled in a separate document 
which is presented as Volume 2. It contains all the assumptions and inter- 
mediate steps which are an important part of the theoretical development but 
which are not included in Volume 1. The main purpose of Volume 2 is to 
provide a study guide or reference for those interested in the theoretical 
aspects of the method of matched asymptotic expansions and/or those who 
may wish to modify or extend the results contained in Volume 1 to' fit some 
particular problem. 

Inasmuch as each volume was written as a separate document, there is a 
certain amount of ove,rlap and cross referencing'beitween the two. Thus the 

t 

reader desiring a more detailed discussion of a particular section in 
Volume 1 need only refer to the corresponding section in Volume 2 and need 
not read through the entire theoretical analysis. 
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ABSTRACT 


Previously published asymptotic solutions for lunar and 
interplanetary trajectories have been modified and combined to 
formulate a general analytical solution to the problem on 
N-bodies. The earlier first-order solutions, derived by the 
method of matched asymptotic expansions, have been extended 
to second order for the purpose of obtaining increased accuracy. 
The derivation of the second-order solution is summarized by 
showing the essential steps, some in functional form. 

The general asymptotic solution has been used as a basis for 
formulating a number of analytical two-point boundary value 
solutions. These include Earth- to-moon, one- and two-impulse 
moon-to-Earth, and interplanetary solutions. Each is presented 
as an explicit analytical solution which does not require interative 
steps to satisfy the boundary conditions. All required formulas 
are presented for each solution. 

Comparisons between the asymptotic solutions and numerical 
integration are shown for several applications. The results show 
that the accuracies of the asymptotic solutions range from an 
order of magnitude better than conic approximations to that of 
numerical integration itself. Also, since no iterations are 
required, the asymptotic boundary value solutions are obtained 
in a fraction of the time required for comparable numerically 
integrated solutions. 

The subject of minimizing the second-order error is discussed , 
and recommendations made for further work directed toward 
achieving a uniform accuracy in all applications. 
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Section 1 
INTRODUCTION 


A number of approximation techniques have recently been proposed for 
calculating N-body trajectories (where N is greater than two). These tech- 
niques include the matched asymptotic expansion (References 1 and 2), 
hybrid patched conic (Reference 3), overlapped conic (Reference 4), multi- 
conic (Reference 5), virtual mass (Reference 6), slowly varying functions 
(Reference 7), and Chebyeshev series (Reference 8). All these techniques 
are claimed to be much faster than numerical integration and considerably 
more accurate than the well known patched-conic approximation. Of all these 
techniques, the matched asymptotic expansion is somewhat unique since it 
represents an analytical solution to the problem of N bodies rather than just 
a numerical scheme for rapid calculation. The analytical nature is useful 
in solving two-point or mixed boundary value problems since, in most 
instances, the solution can be obtained explicity and does not require itera- 
tive steps. 

The N-body problem is one of determining the motion of a body of negligible 
mass subject to the gravitational forces of one primary body and N-2 second- 
ary bodies. The motion of the secondary bodies relative to the primary body 
is assumed to be known. In general, the dominant force on the negligible 
mass body is that of the primary body. However, during a close approach 
of any one of the secondary bodies there is a change in the ordering of the 
dominant and perturbing forces and as a result the problem falls into a class 
known as singular perturbation problems (Reference 9). An approximate 
analytical solution can then be obtained by the method of matched asymptotic . 
expansions. 

Numerical schemes give solutions to this type of problem but they require a 
prescribed state vector at some time t = t in order to uniquely define the 
trajectory. t In many boundary value problems the initial state vector is not 
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known a priori but is only partially prescribed along with some terminal 
conditions. The from er'i cal schemes then require an iterative procedure to 
obtain the unknown part of the initial state vector, i. e. , to solve the two- 
point boundary value problem. 

The asymptotic solution can be formulated to solve the two-point boundary 
value problem directly, i.e. , the unknown part of the initial state vector can 
be obtained without iterations. The solution is formulated as a set of analy- 
tical expressions in the form of asymptotic expansions. Evaluating the 
expressions in a certain sequence gives all the unknown parameters as func- 
tions of the prescribed boundary conditions. The goal of this study was to 
formulate a general, second-order asymptotic solution to the problem of 
N-bodies and to construct from this solution several two -point boundary value 
solutions. This goal can be divided into three specific objectives. 

The first objective was to extend the previously published first-order solu- 
tions to second order. The results of this effort are summarized in 
Section 2, where the N-body differential equation of motion is used as a 
starting point. Section 2 covers the development of the outer and inner solu- 
tions, the overlap domain, the matching, and the fundamental solution. The 
latter gives the relationships between the constants of motion of the outer ! 
solution, where the primary body is dominant, and the inner solution, where 
one of the secondary bodies is dominant. 

The fundamental solution was used to achieve the second goal of this study, 
the formulation of several different asymptotic two-point boundary value 
solutions. These solutions, which can be applied to certain classes of Earth- 
to-moon, moon-to -Earth, and interplanetary trajectories, are presented in 
Section 3. Two versions are presented for each solution, one linear, the 
other nonlinear. In every case at least one of the two solutions satisfies the 
boundary conditions exactly without iterations giving an explicit boundary 
value solution. 

The third objective of the study was to compare the asymptotic boundary 
value solutions with numerical integration. Comparisons for Earth-to-moon, 
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two-impulse moon-to -Earth, interplanetary midcourse, and interplanetary 
trajectories are presented in Section 4. These results show that (1) the 
interplanetary solutions are more accurate than the lunar solutions, (2) mid- 
course solutions are more accurate than those which originate close to one 
body and terminate close to another, (3) the second-order solutions improve 
the first order in some but not all applications, and (4) the computation times 
for the asymptotic solutions are 6 to 150 times faster than for numerical 
integration. 

A discussion of the conclusions obtained from this study and recommenda- 
tions for. further studies are contained in Section 5. 

This study has focused on the application of the method of matched asymptotic 
expansions, and it has assumed that the reader has a certain degree of 
familiarity with the theoretical background. An excellent discussion of the 
basic theory can be found in Reference 9. 

The notation used in this report is a combination of that of Lancaster 
(Reference 1) and Carlson (Reference 2). In general, each parameter is 
defined as it is introduced, but some which have only mathematical meaning 
and serve an intermediary role are defined only by an equation. Scalars 
are written as x or X and vectors as x or _X.. A matrix G(x) and a tensor 
H(x) are also used. In addition, a bar over a parameter indicates that it 
applies specifically to an inner solution. Finally, the order of a particular 
term in an expansion is given by the exponent of the parameter p which 
precedes the term, i. e. , p n is order n or 0(n). 
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Section 2 

SUMMARY OF THE ASYMPTOTIC N-BODY SOLUTION 

The problem of N-bodies for which an asymptotic solution is desired consists 
of finding the motion of a body of negligible mass (hereinafter referred to as 
the particle) under the influence of one primary body and N-2 perturbing 
bodies whose motions relative to the primary body are known. This problem 
requires the solution of the differential equation 


r = _f(r) + F(r, £i ) 


( 2 - 1 ) 


where _r is the position of the particle with respect to the primary body and 

th 

is the position of the i — perturbing body with respect to the primary body. 
The functions f_ and F are defined by 


f ( r) = -r/r' 


( 2 - 2 ) 


N-2 


ZK £, £.) = J Hi [l(£-£i) +1^)] 


(2-3) 


i= 1 


Equations (2-1) to (2-3) are dimensionless; the unit of length is the semi- 
major axis of the orbit of the j— body,* and the unit of time is the period of 
the j — body divided by 2 tt. In dimensionless units then the mass of the pri- 
mary body is unity and of the i— body is |L, which is assumed to be much 
less than one. The origin of the coordinate system is the primary body 
rather than the center of mass and this gives rise to the last term in (2-3). 


*The j — body will be termed the reference body. 
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For lunar trajectories the primary body is the Earth, and the two perturbing 
bodies of interest are the moon and the sun. (Although the p^ for the sun is 
not small compared to the unit mass of the Earth, its contribution is small 
due to the great distance of the sun from the Earth. ) For interplanetary 
trajectories the primary body is the sun and the perturbing bodies the 
planets. 

As long as the particle is not close to one of the perturbing bodies, the 
function F^ in (2-1) is small compared to the other two terms. However, if 
a close approach is made to one of the perturbing bodies, then F becomes 
the dominant force and the problem falls into a class known as singular 
perturbation problems (Reference 9). An approximate analytical solution 
can then be obtained by the method of matched asymptotic expansions. 

The asymptotic solution is formulated by considering two limits of (2-1) and 
then matching the corresponding solutions in an overlap domain. The result 
is termed the fundamental solution and is used to formulate the boundary 
value solutions in Section 3. 


2. 1 OUTER LIMIT 

The outer limit is defined as the limit where r-jx = 0(1) for all i. Then F 
is always small in (2-1), and the solution is assumed to be given by the 
asymptotic expansion 

£(t) = r^t) + H L £ 1 ( t ) + N^£. 2 (t) + 0(p 3 ) (2. 1-1) 


Where the reference mass |i is equal to p^, the dimensionless mass of the 
reference body. Substitution of (2. 1-1) into (2-1) and equating powers of p 
leads to the differential equations for r^, _r and r^. They are 


£i = + Zido'Pi) (2.1-3) 

'±Z = G(r o )r 2 + F^, r y p.) (2. 1-4) 
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where 


N-2 


£1 ( io>Pi ) = 2 M I [i(£ 0 -Pi ) + i ( -Pi ) ] 
i = l 


(2. 1-5) 


N-2 


r z (r o ,r 1 , E .)=iH(r o )r 1 2 + £ (2.1-6) 


i=l 


and 


M. = p./p 


(2. 1-7) 


The function G (x) is a matrix defined by 


3x.x. 

a. = 

ij 5 3 


x x 


( 2 . 1 - 8 ) 


where 6 „ is the Kronecker delta. The function H(x) is a tensor defined by 


15x.x.x, - 

= V — + — C (x. 6., + x. 6 ji_. + x. 6,.) 


-ijk 


■ -j - ik • ~k u ij' 


(2. 1-9 


G and H represent the first and second derivatives in the Taylor series 
expansion of J[( r) about the nominal value _r = r . 


The solutions of (2. 1-2) through (2. 1-4) depend on the initial condition on r_ 
and the corresponding velocity v. These initial conditions can be stated as 


r(t Q ) = £ 0 ( t 0 ) + ^-l^V + H- 2 £ 2 ( t 0 ) (2 - 1_10 ) 

v(t Q ) = v o (t Q ) + pv 1 (t Q ) + p 2 v 2 (t Q ) (2. 1-11) 
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Then the solutions are 


Io (t) = f o {t) ^o (t o ) + 


( 2 . 1 - 12 ) 


jl^t) = A(t, t Q )r 1 (t o ) + B(t 


•‘o ) 2.1 <t o ) + / 


B'(t, t) Fj|T)dT (2, 1-13) 


r 2 (t) = A(.t, t Q )r 2 (t & ) +. B(t, t o )v 2 (t Q ) + J B(t, t ) _F 2 (t )dr 

t 


l <>’ './ 


(2. 1-14) 


The solution for is the standard two-body ellipse resulting from the two- . 
body differential equation (2. 1-2). : The functions f and g Q are infinite 
series in time or can be written in closed form using eccentric anomaly as : 
the independent variable. They are defined in any standard astrodynamics 
textbook and in References 2 and 10.* The solutions for r^ and r_ 2 are made 
up of a homogeneous solution which is simply the propagation of initial devia- 
tions along the two-body solution jr Q , and a particular integral which intro- 
duces the perturbations from two-body motion. The functions A(t, t Q ) and 
B(t, t ) are partial derivative matrices which arise by partitioning the state 
transition matrix 


f'9r Q (t) 3r o (t) 


$ (t. t 


9r (t ) 9v (t ) 
-o o -o o 


3v o (t) 9v 0 (t) 


A(t, t ) B(t, t ) 
o o 


(2. 1-15) 


?r (t ) 9v (t ) 
1-0 o -o o ‘ 


C(t, t Q ) D(t,t o ) 


^-Reference 10 has been included as a second volume of this report. It is 
hereinafter referred to only as Volume 2.. - 
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Many expressions have been derived for the partial derivative matrices in 
various coordinate systems. Some are discussed in Reference 2 and in 
Volume 2, 

The solutions given by (2. 1-12) through(2. 1-14) can be substituted back into 

(2. 1-1) giving a second-order solution for the position £ of the particle. The 

solution is a function only of the initial conditions (2. 1-10) and (2. 1-11) and 

the time histories of the perturbing bodies p.. The initial conditions can be 



chosen so that at some time t = t^ the trajectory passes close to the k— body. 
This introduces another limit. 

The outer solution and its behavior as t approached t^ are shown in 
Figure 1. 

2. 2 . INNER LIMIT 

The inner limit is defined as the limit where £-pv = 0(p ) for some k. This 

^ ^ til 

limit arises when the particle makes a close approach to the k— body. This 
limit requires the change of variables 

" 7 “ " 7 " " * " " : " _ _ CR17 



Figure 1. Outer Solution, Inner Solution, and Overlap Domain 
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( 2 . 2 - 1 ) 


Rk = (£-P k )/n k 


S k = (t_t pk ) ^k 


( 2 . 2 - 2 ) 


where p, is the position of the k — body and t , is the time of pericenter 

— * th P k 

passage of the trajectory about the k — body. The latter can be written 


V = V + ^k T k 


(2. 2-3) 


•and t, chosen as the time at which the two-body outer solution r passes 

k th ^ 

through the center of the k — body, i. e. , at t=t k 


£o (t k> = 2 k (t k ) 


(2. 2-4) 


Substitution of the inner variables (2. 2-1) and (2. 2-2) into (2-1) gives the 
inner differential equation 


d ^k 


= I<*k> + ^-Bk^i) 


(2. 2-5) 


where P is a function defined in Volume 2. If the expansion 


**< S k> = \ 0 < S k> + + ^k 2R k2< S k» + 0 ^k 3 > < 2 - 2 - 6 > 


is substituted into (2. 2-5) and powers of equated the differential equations 


d ^kc 




(2. 2-7) 


d ^kl 


( 2 . 2 - 8 ) 
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1 


d2 ^kZ 

-4- = G ^o^2 + < 2 * 2 ‘ 9 > 

db, 

k 


The solutions of (2. 2-7) through (2. 2-9) depend on the initial conditions on 
and the corresponding velocity V^. These initial conditions can be stated 


as 


iy S ko> = *ko< S ko> 


( 2 . 2 - 10 ) 


^k (S ko> = X ko (S ko ) 


( 2 . 2 - 11 ) 


These initial conditions assume that the perturbations vanish at S, = S, and 

k ko 

that the full solution (2. 2-6) can be represented by at this point. As a 

result the solutions are 


iW S k> = W*W S ko> + *o< S k ) ^o< S ko ) 


(2. 2-12) 


-^fcl^k* 0 


^2< S k> = 



B( s k ^ )G (p k (t k ))R ko u) d e 


(2. 2-13) 


(2. 2-14) 


The solution for R^ 0 is the standard two-body hyperbola resulting from the 
two-body differential equation (2. 2-7). The functions f and g Q are defined 
in many textbooks and in Volume 2. Like their counterparts in the elliptical 
solution of the previous section, they can be written in closed form using 
hyperbolic eccentric anomaly as the independent variable. The function 
B(S , £) is the partial derivative matrix 

1C 


B(S k , | ) 


9 iW s k> 

8 ^ko<6> 


(2. 2-15) 
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The solutions given by (2.2-12) through (2.2-14) can be substituted back into 
(2. 2-6) giving a second-order solution for R^. The solution is shown in 

Figure 1. 


2. 3 OVERLAP DOMAIN 

The outer and inner solutions are functions of the vector constants given by 
(2. 1-11), (2. 1-12), (2. 2-10), and (2. 2-11). For a trajectory which is con- 
tinuous from the domain of the outer solution to the domain of the inner 
solution, these constants are not all independent. In a region designated as 
the overlap domain, the outer and inner solutions must exhibit certain 
Similarities, i. e. , both solutions must represent the trajectory in this 
domain. This characteristic makes it possible to determine explicit rela- 
tionships between the constants of the two solutions. A representation of the 
overlap domain is shown in Figure 1. 


The overlap domain is defined as the domain of the intermediate limit, 
i. e. , the domain where t-t ^ = (Mp 01 ) with 0<a<l. This limit is formally 
defined by introducing the intermediate variable 


(T, = (t-t . )/p!* 0<a <a:«x: <1 (2. 3-1) 

K p.K ^ O I 

•i 

If a = 0 then (2, 3-1) simply shifts the time scale of the outer solution to a 

new origin. If a = 1, giving the inner time. Within the range 

sn<rn-,< cr , is then intermediate to the outer and inner times. The values 
°o^“^“l k 

of o q and must be determined from the matching. 


The outer solution is a function of t and replacing t 
and solving for t-t^ gives 


pk 


in (2. 3-1) by (2. 2-3) 


* _t k " k + H k T k (2.3-2 

Since p^_ is small, (2, 3-2) indicates that the outer solution must be expanded 
to about t = t^ in order to determine its behavior in the overlap domain. 

This expansion is derived in Section All of Volume 2. 
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The outer expansion can be summarized in function form by the following 
expression: 


£ = £(t-t k , 


n; £.(* 0 ). v(t o ); 


-k' : ’ -k*’ J»k*’ — k : ) 


(2. 3-3) 


The position vector _r, when t-t^ is small, is a function of both the initial 

conditions at t and four constants, andp^*, which represent 

the first- and second-order deviations from two-body motion over the 

interval t <t<t, . These constants are discussed in detail in Subsection 3. 8. 
o k 


The inner solution is a function of and comparing (2.3-1) with (2. 2-2) 
yields 




(2. 3-4) 


Since a - 1 < 0, (2. 3-3) indicates that the inner solution must be expanded 
for S^ large. This expansion is derived in Section A 12 of Volume 2. 

The inner expansion can be summarized in functional form by the following 
expression: 


-^k = ^k (S k’ -ook' — k ; — kf) ( 2 ' 3 - 5 > 

The position vector R^, when S^ is large, is a function of X<x>k’ ^yper- 
bolic excess velocity, _L a vector function of the orbital elements, and A^, 
which represents the second-order deviation from two-body motion far out 
on the asymptote of the zeroth-order hyperbola. These constants are dis- 
cussed further in Subsection 3. 8. 



The expansion is obtained by a Taylor series in Section All of Volume 2. 

It may be summarized simply as the function 

P k = Pk (t ' t k ) (2 * 3 ~ 6) 
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2. 4 MATCHING 

For the outer and inner solutions to match, they must be in terms of a 
common independent variable. The expansion of the outer solution near 
t = t^, summarized by (2. 3-3), can be written in terms of and using 
(2. 3-2) giving 

- ~ (r k + M 'k T k’ fi: - ( t 0 ) ’^ (t o ); -k‘ S — k : ’ — k ' ’ -k") = j|_l (2.4-1) 

The position t _ in terms of the inner solution is found from (2. 2-1), i. e. , 


-t = Ek +l *A 


(2.4-2) 


th 


The expansion for the position of the k — body near t = t^, given by (2. 3-6), 
can also be written in terms of cr^ and giving 


Ek = PkK "k + IVk)- *2 


(2.4-3) 


Finally the expansion for when is large, obtained from (2. 3-5), can be. 
written in terms of o- and using (2. 3-3) giving 


^k “ ^k(^k ^k’^k 1 — ®k’— k : — k2) = —3 


(2.4-4) 


Substituting (2.4-3) and (2.4-4) into (2.4-2) gives 

£ = + ^kiL3 { ' 2 " 4 "^ 

Simply stated, the matching requires that the difference between the outer 
solution, as given by (2.4-1), and the inner solution, as given by (2.4-5), 
must be vanishingly small in some appropriate limit. Cole (Reference 9) 
states this limit as 
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lim 


^k ~ ° 


constant 


6 (nJ 


= 0 


(2. 4-6) 


where € (p^) is a gauge function. For a second-order theory e ( h-^.) is most 
easily chosen to be . 


In Section A 14 of Volume 2, it is shown that this limit exists only if 


a = 2/5 


(2.4-7) 


a 1 = 1/2 (2.4-8) 

Thus the overlap domain is a region of order p Q where 2/5 < a < 1/2 and 
a = 1/2 is not included. This is a result of the second-order solution inas- 
much as Carlson (Reference 2) showed that the first-order solution can be 
matched with a = 1/2. 

This was an assumption in his derivation and not a result of applying a 
rigorous matching requirement such as (2.4-6). The present results show 
that his approach to matching will not work for second order, i. e. , certain 
terms which are singular in the limit (2. 4-6) can only be eliminated if 
a < 1/2. 


2. 5 FUNDAMENTAL SOLUTION 

The complete matching process is discussed in Sections A14 to A17, and 
Section B1 of Volume 2. The result, summarized in one six-component 
state vector equation which will be called the fundamental solution, is 


r i(t o ) + H£ 2 (t o ) 

\ V!(t 0 )+pv 2 (t o ) 


X 


L_k 


Xk + ^ik 


* <V V < 


-1 


* (^cok-^kV 


(2. 5-1) 
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where 



(2.5-2) 

(2.5-3) 
(2. 5-4) 
(2. 5-5) 



(2. 5-6) 


and, except in the log term, p has been eliminated in favor of the reference 
mass (jl (which may be equal to if the k — body is also the reference body 

used to non-dimensionalize the differential equations). 


Equation (2. 5-1) is a relation between the constants of the outer and inner 

solutions. The only constants which do not appear explicitly are the initial 

position and velocity of the zeroth-order outer solution, r (t ) and v (t ). 
r 7 -oo-oo 

They must be chosen to make the zeroth-order ellipse intersect the position 

of the k— body at t = t^, i. e. , to satisfy (2, 2-4), They then enter implicitly 

through the relative velocity _V, which is the difference between the zeroth- 

^ th 

order velocity and the velocity of the k — body at t = t^ (and should not be 
confused with the inner* time -dependent velocity V^S^.)), 

Equation (2. 5-1) can be used to solve either initial or boundary value prob- 
lems. The initial value solution is discussed in Section A of Volume 2. 
It is the boundary value solution which is of interest in this study, and the 
applications of (2. 5-1) are discussed in the following subsections. 
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2. 6 COMMENTS 

The asymptotic solution presented here is similar, when second-order terms 
are ignored, to the first-order solution derived by Carlson (Reference 2). 
The obvious differences between the two first-order solutions are (1) the use 
of dimensionless variables, (2) isolating the small parameter p so that it 
appears explicitly, and (3) using the vector _L^ as one of the constants of the 

v 

inner solution rather than the standard impact parameter vector. 

The two solutions are numerically equivalent when applied to an initial value 
problem. However, using the vector does result in a mathematically 

different solution when the fundamental solution is applied to boundary value 
problems. This is because the use of the impact parameter vector results 
in boundary value solutions which satisfy the boundary conditions in a "best" 
sense while the vector results in solutions satisfying the boundary con- 
ditions exactly. The two vectors are related by (cf. Volume 2) 


4 - 


^k + (t V s k>^.k 


( 2 . 6 - 1 ) 


where and n^ are defined in Subsection 3. 8. 


Next, it should be noted that the second-order terms add considerable com- 
plexity to the solution although such complexity is not apparent in this 
section. Some of the complexity can be seen from the formulas in Sub- 
section 3. 8, but it is necessary to follow the derivation in Volume 2 to really 
appreciate just how much complexity is actually added. The amount of 
algebra necessary to extend the solution to a higher order would probably be 
prohibitive and the result somewhat unmanageable. The first-order solution 
contains the 3x3 gravity gradient matrix G, while the second-order" solution 
contains the 3x3x3 tensor H. Each succeeding order adds a tensor of higher 

order. If the dimensions of the tensor are used as a measure of the com- 

th 

plexity of the solution, then an n — order solution has a complexity of order 
3 n+l , 
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Finally, it should be pointed out that the form of the fundamental solution is 
not unique. The matching results contained in the fundamental solution are 
similar, but not identical, to those of Carlson (Reference 2). Differences 
which are not immediately obvious are due to the fact that an asymptotic 
expansion of a given function is not unique. Other expansions can be formu- 
lated to represent the same function but actually appear as different expan- 
sions. When the individual expressions which result from the matching are 
combined to form the fundamental solution, there are several ways in which 

such a combination can occur. Thus for a second-order solution the error 

3 ; 

in each case may be. order p , but the actual value of the error may differ, 

3 3 

i. e. , for one case it may be 3p while for another case it may be 0. 5|jl . 

This aspect is discussed further in Subsection 4. 6. 


Section 3 


ASYMPTOTIC BOUNDARY VALUE SOLUTIONS 

The boundary value solutions presented in this section are of three general 
types: (1) trajectories which originate at some known position relative to 

the primary body and terminate at a fixed pericenter radius, inclination, and 
time at one of the perturbing bodies {cf. , Subsections 3. 2, 3. 3, and 3.6), 

(2) trajectories which originate at some known position close to one of the 
perturbing bodies and terminate close to the primary body with fixed entry 
conditions (cf. , Subsections 3.4 and 3. 5), and (3) trajectories which originate 
close to one perturbing body and terminate close to another with fixed 
pericenter radius, inclination, and time at each end (cf. , Subsection 3. 7). 

Each of the boundary value solutions evolves from (2. 5-1),* and each requires 
at least one solution of a Lambert problem to establish the zeroth-order 
outer solution. The two types of Lambert solutions which are required are 
discussed in Subsection 3. 1. Subsections 3. 2 through 3,7 present the 
various boundary value solutions, and finally Subsection 3. 8 gives formulas 
for evaluating all the constants which appear in the boundary value solutions. 

The sections which follow contain only the end results of the boundary value 
solutions. More detailed discussions and the steps necessary to go from 
(2. 5-1) to each solution are contained in Volume 2. 

3. 1 THE LAMBERT PROBLEM 

The standard Lambert problem is one of finding the two-body solution which 
connects two known position vectors in a fixed time of flight. If the two 


^Except for the two-impulse moon-to-Earth solution for which two different 
types of solutions have been derived, one of which does not evolve from 
(2. 5-1). It is this latter solution which is presented in Subsection 3. 5. 
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position vectors are and x.^ and the time of flight t f , then Lambert's 
theorem states that 


t^ - t^(a, + x^, c) - t^ - tj (3. 1-1) 

where a is the semimajor axis which is unknown and c is the chord length 
between x^ and x^- The chord length is found from the law of cosines, i. e. , x 


2 

c 


x. 


+ x^ + 2x^ x^ cos 0 


'12 


(3.1-2) 


where 9^ is the central angle between x^ and x^. 


T 


An iterative solution is required to determine the semimajor axis. Once it 
is known, the velocities x^ and x^ and the solution x(t) can be obtained. Many 
techniques have been proposed for solving the Lambert problem. One such •> 
'method is discussed by Battin (Reference 11). 


The standard Lambert problem requires that the two vectors and x^ be 
given. The zeroth-order solutions used in the Earth-to-moon and inter- 
planetary solutions are of this type. The zeroth-order moon-to-Earth 
solutions however do not rely on a given position vector at the Earth. 

Instead, entry conditions of radius and flight path angle are prescribed at a 
given time. In addition, the trajectory is to satisfy a prescribed inclination. 
The solution for the semimajor axis is now more difficult, since 0 ^ i n 
(3. 1-2) is not known a priori. This angle is the difference between the true 
anomalies at the endpoints, i. e. , 




(3. 1-3) 


where 


- 1 

cos 


a (1 - e 


ex. 



(3. 1-4) 
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The eccentricity e is a function of a, x^ and the flight path angle "Y^, 
measured from the local horizontal. The relationship is 

/a (3.1-6) 



The modified Lambert problem requires the simultaneous solution of 
(3. 1-1) through (3. 1-6) for a, c, 0^, ^ anc ^ e " ® nce these parameters 

are determined they, along with the prescribed inclination, are sufficient to 
solve for the velocities and and the time-dependent solution x(t). The 
solution is discussed in detail in Section B4. 1 of Volume 2, and a similar 
problem is discussed in Reference 11. 

3.2 EARTH- TO-MOON SOLUTION 

The Earth-to-moon problem is one in which the target body is the moon. 

The moon should also be the reference body, therefore 

k = M (3.2-1) 

M- = (3. 2-2) 

The simplest Earth-to-moon boundary value problem is shown in Figure 2. 

The initial time, t , the initial position relative to the earth, r(t ), and the 

o _ — o 

pericynthion radius, p^, inclination, i and time, tpj^, are all pre- 
scribed. The initial velocity relative to the earth, v(t Q ), is unknown and 
must be determined from the fundamental solution. In order to evaluate the 
solution, an ephemeris is required giving the position and velocity of the 
moon and the position of the sun in Cartesian coordinates with origin at the 
Earth. 
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Figure 2. Earth-to-Moon Solution 


From Figure 2 it can be seen that the zeroth-order ellipse, r_ Q (t), coincides 
with the higher order solution, jr(t), at t = t Q . Therefore, in (2. 1-10) let 


“ -2 {t o ] = 0 


Then 


r (t ) = r(t_) 

o o — o 


defining the initial position of the zeroth-order ellipse.'' 
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From (2. 2-3) 


t M ~ ‘PM ^ T M 


(3. 2-5) 


where is arbitrary and can be put equal to zero without loss of generality. 
Non- zero values of t^. simply causae a change in t^ when is held con- 

stant. The final'positiori of the zeroth-orde'r ellipse comes from (2. 2-4) 




(3. 2-6) 


where p^ is the position of the moon obtained from the ephemeris. The two 

position vectors, _r (t Q ) and _r Q (tj^), define a standard Lambert problem of 

the type discussed in Subsection 3. 1. Solution of the problem gives r^(t), 

shown as the dashed line in Figure 2, and the initial and final zeroth-order 

velocities, v (t ) and v (b ,). The latter is used to define the relative 
— o o — o M 

velocity 


M " -o (t M ) 


(3. 2-7) 


where is the velocity of the moon. 


Now let the initial velocity perturbation be 


6v(t ) = v,(t ) + p v_(t ) 

— O — i O — c. O 


(3. 2-8) 


Then the initial velocity at t is 

o 


v(t ) = v (t ) + p .6 v(t ) 

— o' — o o — o 


(3. 2-9) 


while the excess velocity at the moon can be written 


V ™ = V + p 6 v 

— oo M — M — 


00 


(3. 2-10) 
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The perturbation terms in (3. 2-9) and (3. 2-10) are obtained from the 
fundamental solution. They are 


6v(t o ) = 




+ -M + ^ -M* 


(3. 2-11) 


6V, = D(t M , t o ) 6 v(t o ) - 6 m -p n M ‘ (3.2-12) 

Equations (3. 2-9) through (3. 2-12) constitute a linear solution to the 
boundary value problem. 

O' 

Since enters (3.2-11) through and since 6 v(t ) appears in (3.2-12) 

they are not explicit relations but must be solved in a sequence using the 
zeroth-, first-, and second-order terms successively. The zeroth-order 
approximation is obtained by putting p = 0 in (3. 2-9) and (3. 2-10). The 
first-order approximation is obtained by putting p = 0 in (3. 2-11) and 
(3. 2-12) and using the zeroth-order value of V ^ to evaluate (3. 2-11). The 
second-order approximation is obtained by evaluating (3. 2-11) and (3. 2-12) 
with p 4 0 and using the first-order V ^ in (3. 2-11). 

Combining (3. 2-9) with the prescribed value of £(t Q ) gives a complete set 
of initial conditions for a trajectory satisfying all the conditions of the 
boundary value problem. Combining (3. 2-10) with the prescribed values of 
pericynthion radius and inclination gives a complete set of terminal condi- 
tions as shown in Subsection 3. 8. 

An alternate solution, called the nonlinear solution (Reference 2), can be 
obtained from the solutions of a sequence of Lambert problems defined by 

C 

the position vectors 

r'(t ) = r (t ) (3.2-13) 

- — o o — o o . 


4 (t M> • + 


(3.2-14) 
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where 


6 i'W = Sm+Im + i'Im (3 - 2 - 15 > 

Solution of the Lambert problems gives the initial and final velocities, 

v 1 (t ) and v 1 (t. Then the initial velocity replacing (3. 2-9) is 
— o o — o M 

V(t j = v'(t ) (3. 2-16) 

— o o o 

; r 

and the excess velocity replacing (3. 2-10) is 

— ocM = ^M + ^lo (3.2-17) 

where 

Ym = - £m<W < 3 - Z - 18 » 

‘^1 = - ® M - I 3 - 2 - 1 ’) 

Again the solution requires a sequence- of steps. The zeroth-order 

approximation is obtained by putting p = 0 is (3. 2-14) and (3. 2-17) and is 

identical to the zeroth-order linear solution. The first-order approximation 

is obtained by putting p = 0 in (3. 2-15) and (3. 2-19) and using the zeroth- 

order V (3. 2-15). And the second-order approximation is obtained by 

using the first-order V ^ , in (3. 2-15). The nonlinear solution is shown in 

— ooM. 

Figure 3. 

The first- and second-order nonlinear solutions will be slightly different 
from their linear counterparts since they include nonlinear effects in the 
zeroth-order solution which are not contained in the B and D partial deriv- 
ative matrices used in the linear solutions. 


25 


The constants "Y^, an< ^ are fixed through each step of both the 

linear and nonlinear solutions. The function however, depends on 

V and must be evaluated for each of the zeroth- , first- ,and second-order 
approximations. Formulas for calculating all of the constants are in 
Subsection 3. 8. 

3. 3 EARTH-TO-MOON MIDCOURSE SOLUTION 

In the previous section, the initial position, _r(t o >, was implicitly assumed to 

be close to the Earth. The same analysis may also be used for a midcourse 

maneuver where the position, _r(t ), represents a point between the Earth and 

the moon, as shown in Figure 4. The velocity just prior to the midcourse 

+ 

maneuver is v(t ) and after the maneuver it is v(t ). Therefore, the 
— o — o 

midcourse velocity correction is 

Av(t o ) = V(t+) - v(t;> 

= v (t + ) - v(t ) + p6v(t + ) (3. 3-1) 

— o o — o 1 — o 
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The solution of Subsection 3. 2 can be used to calculate v (t + ) and 6v(t + ) 

— o o — o 

and since _r(t Q ) and v(t Q ) are known, (3. 3-1) gives an analytical expression 
for the midcourse velocity correction. 


3.4 ONE-IMPULSE MOON- TO-EARTH SOLUTION 

In the moon-to-Earth problem, the moon becomes the launch body and is 
also the reference body. Therefore, as in Subsection 3. 2, 

k = M (3.4-1) 

p = P M (3.4-2) 

The boundary value problem is shown in Figure 5. The initial time, t^, 

the initial position relative to the moon, an< 3 the entry time, t , 

radius, r , flight path angle, y , and inclination, i , are all prescribed, 
e £ e 

The initial velocity relative to the moon, is unknown and must be 
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determined using the fundamental solution. An ephemeris like that used for 
the Earth-to-moon solution is required. 

Since entry conditions rather than a fixed position vector are prescribed at 
Earth, this solution requires the solution of the modified Lambert problem 
discussed in Section 3. 1. The initial position of the zeroth-order ellipse 
comes from (2. 2-4) 


io (t M> = PM^M 1 


where p M is the position of the moon obtained from the ephemeris and, 
for convenience, 



(3.4-4) 
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The final position r (t ) must be determined from the solution of the 
modified Lambert problem along with £^(t), shown as the dashed line in 
Figure 5, and the initial and final zeroth-order velocities, V Q (tj^) and 
v (t ). The initial zeroth-order velocity is used to define the relative 
velocity 

— M = W -£m ( V (3 ' 4 - 5) 

where p^. is the velocity of the moon. 

The zeroth-order solution satisfies the entry conditons exactly, therefore 
any perturbations at t = t will cause the trajectory to deviate slightly from 
the prescribed conditions. As shown in Figure 5, the position perturba- 
tion can be made to vanish, i. e. , 

r^) = £ 2 (t o ) = 0 (3.4-6) 

so that 


r(tj = r (t ) 
— o — o o 


(3.4-7) 


The velocity perturbation cannot vanish without overly constraining the 
problem, therefore, 


6 v(tg) = v^tg) + p v 2 (t e ) i 0 


(3.4-8) 


so that entry velocity is 


= X,( f e) + l x6 l( t e ) 


(3.4-9) 


The hyperbolic excess velocity is given by 


— coM = + 


(3. 4-10) 
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The perturbation terms in (3.4-9) and (3.4-10) are obtained from the 
fundamental solution. They are 


6v(t e ) = B(t M , t e ) 

( -M + -M + K -^M ) 

(3.4-11) 

= D(t M , y 

6v < t e) - k M - n M 

(3. 4-12) 


Using the initial position R w , and the excess velocity V a new para- 

— Ml — oo M 

meter X-^ is defined by 


X M = 1+4 


r mi V ooM 


1 + cos 


I R V ' 

-Ml • — ooM 

V R M1 V ooM , 


-1 


(3. 4-13) 


Then the initial velocity at the moon is 



(3.4-14) 


Equations (3.4-9) through (3.4-14) constitute a linear solution to the 
boundary value problem. They must be solved in a certain sequence to 
obtain zeroth-, first-, and second-order approximations (cf. Section 3. 2 for 
a discussion of the steps involved). Equations (3.4-13) and (3.4-14) must be 
included in the sequence since the constants in (3.4-11) and (3.4-12) are 

functions of R,., and V..,,. 

—Ml —Ml 

Combining (3.4-14) with the prescribed value of Rj^ gives a complete set of 
initial conditions (note that they are inner variables and that to obtain dimen- 
sional values Rj^ must be multiplied by p as well as the appropriate dimen- 
sional length scale) for a trajectory satisfying, to zeroth order,, all the 
conditions of the boundary value problem. The terminal state is given by 
(3.4-7) and (3.4-9). The prescribed entry radius is satisfied exactly 
because of (3.4-6), while the flight path angle and inclination differ by 
order p from the prescribed values due to (3.4-8). 
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An alternative, nonlinear solution can also be obtained in a manner similar 
to that discussed in Subsection 3. 2 for Earth-to-moon trajectories. ■ The 
new Lambert problem is defined by the position vectors 


where 


r o (t M } = + ^ 6 ^< t M ) 


(3. 4-15) 


r (t ) 
— o e 



(3. 4-16) 


6 -E<V = + Im + >* i M (3 - 4 - 17 > 

Equation (3.4-15) defines the fixed position vector for the modified Lambert 

problem and (3.4-16) represents the prescribed entry conditions. Solution 

of the new Lambert problem gives the final position, _r ' (t ), and the initial 

i i o e 

and final velocities, v (L .) and v (t ). The excess velocity is 

— o M — o o 

— ooM = + < 3 - 4 ' 18 > 

where 

' %. = . ' (3.4-19) 

— CO = ^ — M (3.4-20) 

The initial velocity relative to the moon is still defined by (3. 4-14). This 
solution requires the same sequence of steps as the nonlinear Earth-to-moon 
solution discussed in Subsection 3. 2 except that the Lambert problems are of 
the modified rather than standard type. 


For each of the zeroth-, first-, and second-order approximations the entry 
velocity is given by 


xiy = 


(3. 4-21) 
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i. e. , the entry velocity is the velocity of the modified Lambert solution. - 
Since all of the Lambert solutions satisfy the entry constraints, the nonlinear 
solution satisfies these constraints to any order. Thus the nonlinear solu- 
tion has the advantage that it satisfies the entry boundary conditions exactly 
rather than to zeroth order as the linear solution does. The nonlinear solu- 
tion is shown in Figure 6. 

The constants Y_ 6-^. and rj_ are again fixed through each step. The 

function ,, however, depends on R. and V , . . and must be evaluated for 
each approximation. In addition, because of (3.4-4), t / 0 and must be 
calculated with the other constants. This is expected since will not, in 

general, be normal to R^^, i* e., the initial position is not pericenter. The 
parameter is a measure of the time between t^ and the time of pericenter 
passage. Formulas for calculating all of the constants are- given in 
Subsection 3. 8. , 


The initial velocity is the velocity after the impulse. The notation 

represents the velocity, at t 4 , i. e., 


— Ml = -Vl) ( 3 - 4 - 22 > 

If V^.(t.) is the velocity before the impulse then the single impulse is 
given by 

AV = v M (t+) - v M (t-) 

a ^ M 1 ’^ M (t'l) < 3 - 4 - 23 > 

3. 5 TWO-IMPULSE MOON-TO-EARTH SOLUTION 

The two-impulse moon-to-Earth problem is similar to the one-impulse 
problem discussed in the previous section except that the initial velocity is 
assumed to be known and does not give a trajectory satisfying the prescribed 
entry conditions. A second impulse is applied at some time prior to reach- 
ing the moon's sphere of influence resulting in a trans-Earth trajectory 
which does satisfy the entry conditions. 

The boundary value problem is shown in Figure 7. The initial time, tj, the 
initial position and velocity relative to the moon, Rj^j and Vj^j, an initial 
impulse, Ij, along the current velocity vector at tj, the time of the second 
impulse, t2, and the entry time, t e , radius, r e , flight path angle, "V e , and 
inclination, i e , are all prescribed. The magnitude and direction of the 
second impulse are unknown and must be determined from the asymptotic 
solution. 

The solution may or may not be derived from the fundamental solution. If the 
second impulse occurs well outside the moon's sphere of influence then the 
fundamental solution, derived from the matching process, should be used for 
the part of the solution between the first and second impulses. However, 
when the second impulse occurs near or inside the sphere of influence, 
Carlson (Reference 2) has shown numerically that a perturbed hyperbola is 
more accurate than the fundamental solution. His results have been verified 
theoretically by considering the order of magnitudes of all the error. 
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Figure 7. Two-Impulse Moon-to-Earth Solution 


terms at various points out to and slightly beyond the sphere of influence 

( 

(cf. Section B4. 3 of Volume 2). 

The solution described here utilizes a perturbed hyperbola between the 
first and second impulses and a perturbed ellipse between the second 
impulse and Earth entry. The solution therefore is not a matched solution 
but simply two asymptotic solutions joined at the point of the second impulse. 

The initial position and velocity (in inner variables) are 


-M^l* _ -Ml 


" -Ml 
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where' indicates the instant just prior to the first impulse which is 
defined as 


(3. 5-3) 

The position after the first impulse is still given by (3. 5-1), but the velocity 
is now , 

i 

' V M (S;> = (1+ Ij) V M (S-j) (3.5-4) 

i. e. , the post-impulse velocity is parallel to the initial velocity. The 
impulse 1^ must be chosen to make 


V .M = V M< S t>- 2 /R m' S 1»*° 


(3. 5-5) 


i. e. , 


, Ii2 (V2/R m ( s i))/V m (S-) - 1 


(3. 5-6) 


This is needed to guarantee a hyperbolic orbit after the impulse. 
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The post-impulse position and velocity uniquely define a zeroth-order 
hyperbolic trajectory about the moon. Position and velocity, Rj^ q and 
as a function of the inner time S, can be determined from standard formulas 
such as (2. 2-11). In addition the position and velocity can be used to deter- 
mine the eccentricity and mean motion of the hyperbola as well as the initial 
eccentric anomaly F^ [cf. Carlson (Reference 2)], Volume 2, or 
Subsection 3. 8]. Then the initial inner time is 

Sj = (e sinh Fj - F^/n (3.5-7) 


and the inner time of the second impulse is 


S 2 “ S 1 + ^2 t l^^ i 


(3.5-8) 


The position and velocity at including the perturbations due to the 
Earth, are 


*M< S 2> 


= *Mo (S 2 )+ ^ G M 


R> . (S,)^ 
— Mo' 2 2 


s 3 

Xmo |S Z> T 


+ a 3 s* 


(3.5-9) 


V M ( S 2> = - V Mo' S 2^ G 


M 


*Mo( S 2> S 2 " 


-Mo (S 2 ) 2 


3 3 

+ 4p A 3 


(3. 5-10) 


where G m and are defined in Subsection 3. 8. Transforming to Earth- 
centered coordinates using (2.2-1) and (2,2-2) gives 


-^2* ~ -2m ^2* + ^ -M ^ S 2^ 


(3. 5-11) 
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V( t -) = (t 2 ) + v M (S 2 ) 


(3. 5-12) 


The solution between t^ and t is given by (2. 1-12), (2. 1-13) and (2. 1-14) 

with t = t and t = t 0 . 
o e 2 

The zeroth-order solution is found from a modified Lambert problem with the 
initial position given by 


r o (t 2 ) = r(t 2 ) (3.5-13) 

and the terminal position defined by 

r (t ) = r (r , y , i ) (3.5-14) 

-o e -o e e e 

where (3. 5-14) represents, the prescribed entry conditions. Solution of the 
modified Lambert problem gives t q (t & ) , r Q (t) shown as the dashed line in 
Figure 7, and the zeroth-order velocities, v (tt) and v (t ). 

“O l “06 

As shown in Figure 7, the position perturbation at t g can be made to vanish 
giving 


^l (t e } = —2 ^ = ° 


(3. 5-15) 


so that 


r (t ) = r (t ) 
— e' -o ' e 


(3. 5-16) 


The velocity perturbation 

6v(t e ) - Vj (t e ) +p v 2 (t e ) ^ O (3. 5-17) 
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does not vanish, so the entry velocity is 

1 v(t e ) = v q (t e ) + p 6 v (t e ) 


(3. 5-18) 


By combining (3. 5-10) and (3. 5-11) with (2. 1-12) through (2. 1-14), the 
entry velocity perturbation and the velocity impulse at t 2 can be written 

6v(t e ) = -B(t 2 , t e ) _1 [r iQ (t 2 , t e ) + p r 20 (t 2 , t e ) J (3.5-19) 

AV 2 = y o (t 2 ) - v(t 2 ) + jx £r n (t 2> t e ) + D (t 2 , t e ) 6v(t e ) j 

+ P 2 X 21 (t 2 , t e ) (3.5-20) 


where^jg, _T^, r 2Q andJT,^ are cons * an * s defined in Subsection 3.8. Equa- 
tions (3. 5-17) through (3.5-19) constitute a linear solution to the boundary 
value problem. Equation (3. 5-18) gives the deviation from the prescribed 
entry conditions of Y g and i [ r e is satisfied due to (3.5-14)j and (3.5-19) 
gives the velocity impulse which satisfies the entry conditions to zeroth 
order. Note that AV 2 is a function of t 2> of ^Ml’ anc ^i (through 

v(t->)), of r , \ and i (through v (tt)), and the perturbations due to the 
moon and sun (through the T's), 


Therefore AV 2 is a function of all the 


boundary conditions. 


An alternative, nonlinear solution can be obtained by solving a modified 
Lambert problem between the position’vector s 


r o (t 2 ) " £. (* “P £io ^2’ t e> * M- J 20 ^2’ (3.5-21) 


r t ) = r ' (r , , i P ) 

— o e ' — o e *e e 


(3. 5-22) 
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where r'(t~.) is the fixed position vector and r ’ (t ^represents the prescribed 
-o 2 -o e 

boundary conditions. Solution of the new modified Lambert problem gives 
the velocities v '(t±) and v '(t ). The impulsive velocity becomes 

— O l> ~ O 6 

AV 2 = v;(tj) -y(t 2 ) +p E n (t 2 , t e ) +p 2 r 21 (t 2 , t e ) (3.5-23) 

and again, AV 2 is a function of all the boundary conditions. 

The entry velocity is 

v(t ? ) = v^ (t e ) (3.5-24) 

and since v (t ) comes from a modified Lambert solution, the nonlinear 
-o e 

solution satisfies the entry conditions exactly to any order. Therefore, 
just as in the one-impulse case, the nonlinear solution has the advantage 
of satisfying the boundary value problem exactly rather than to zeroth order 
as the linear solution does. The nonlinear solution is shown in Figure 8. 

CR17 



Figure 8. Nonlinear Version of Two-Impulse Moon-to-Earth Solution 
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3.6 INTERPLANETARY MIDCOURSE SOLUTION 

The interplanetary midcourse solution is similar to the Earth-to-moon 
midcourse solution and therefore similar to the Earth-to-moon solution. 

The main difference is that this is an N-body solution where motion takes 
place in heliocentric rather than geocentric space. N may vary from 
three up to eleven, the latter being the entire solar system plus the particle. 
The target body (one of the nine planets) may also be the reference body. 
Letting T indicate the target body gives 

k = T (3.6-1) 


P = (J-r 


(3.6-2) 


The boundary value problem is shown in Figure 9. The initial time, t , the 
initial position relative to the sun, _r(t Q ), and the pericenter radius, p 
inclination, i^,, and time, t^^, at the target planet are all prescribed. 
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The initial velocity, v(t Q ), is unknown. An ephemeris giving positions of 
the N-2 planets in Cartesian coordinates relative to the sun is required to 
obtain the solution. 

The solution is identical to the Earth-to -moon solution described in Sub- 
section 3.2. That is, .the initial velocity v(t ) is found from the equations of 
Subsection 3. 2 with M replaced by T and with p equal to rather than p 

The midcourse velocity correction is found from (3.3-1). Both the linear 
and nonlinear solutions are applicable to this problem although in actual 
practice the linear may be sufficient, as conditions favoring the nonlinear 
solution are not as likely to occur as they are in the Earth-to-moon solution. 

3.7 INTERPLANETARY SOLUTION 

The interplanetary solution is similar to the previous solution in that both 
apply to the general N-body problem. There is a fundamental difference 
between this solution and all the others in that the trajectory passes close 
to two perturbing bodies, one at the launch end of the trajectory and one at 
the target end. Letting L indicate the launch body and T the target body 
gives 

•k ( ' ) “ = L (3.7-1) 

k (2) = T ' (3.7-2) 

The reference body may be either L or T, but to be consistent with the 
previous solution let 

K = P T (3.7-3) 

The interplanetary boundary value problem is shown in Figure 10. The per- 
icenter radius, p^, inclination, i^, and time, t are prescribed at both 
the launch planet, k = L, and the target planet, k = T. The hyperbolic excess 
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velocities relative to the sun as well as the pericenter positions and 
velocities relative to the respective bodies are unknown. The solution 
requires an ephemeris giving planetary positions relative to the sun. 


From (2. 2-3) 

(3.7-4) 

t T = \t “ ^ T T (3. 7-5) 

where, as in the Earth-to-moon solution, the t's are arbitrary and can be 
set equal to zero without loss of generality. 
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From Figure 10 it can be seen that the zeroth-order ellipse, r Q (t), passes 
through the launch planet at t = t^ and through the target planet at t = t^,, 
therefore 

?o (t L> - El/V (3 - 7 ' 6) 

: 0 (t T> = Pt (t T> (3.7-7) 

where p^ and p^, are the positions obtained from the ephemeris. 

The two position vectors, r (t T ) and r (t„), define a Lambert problem and the 

solution gives i- q ( t ) , shown as the dashed line in Figure 10, and the initial 

and final zeroth-order velocities, v (t T ) and v (t™). These velocities are 

- O Li - O 1 

used to define the relative velocities. 

Xl = v D <t L > - p L (‘ L ) 

— T = -o ^T^ ' -Et ^T^ 

where p^ and are velocities obtained from the ephemeris. 

Now let the hyperbolic excess velocities be defined by 

-™L = + 6 ^«>L (3.7-10) 

Vcoj, = Vrp + p. 6Vco-p (3.7-11) 


(3.7-8) 

(3.7-9) 
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The perturbation terms in (3.7-10) and (3.7-11) are obtained from two 
fundamental solutions (Reference 12), one with k = L and one with k = T, 
which can be combined to give 



B (t T , t L ) _1 Jjf T + V T + p C T - A (t T , t L ) (jf L + Y l + p C L )J 
- 6 L - (3.7-12) 


6Veo T 


= C (t T , t L ) tf L + V L + p i L ) + p D (t T , t L ) (6 Voo l + 6. l + p a L ) 


- ~ P 


(3.7-13) 


The position and velocity at pericenter relative to either L or T are given 
by the inner variables 


R 


pk 


— pk - - — k 

\ e k 


i 




1 + e kV /2 3ik. 

9 wk 


(3.7-14) 


(3. 7-15) 


Equations (3. 7- 10) through (3. 7- 1 5) constitute a linear solution. of the 
boundary value problem. Since V^, and Va, enter the right-hand sides of 
(3.7-12) and (3. 7-13) through and the relations are not explicit but 
must be solved in a sequence using the zeroth-, first-, and second-order 
terms successively. The zeroth-order approximation is obtained by putting 
p = 0 in (3.7-10) and (3.7-11). The first-order approximation is obtained by 
putting p = 0 in (3.7-12) and (3.7-13) (but not p L or p T ) and using the zeroth- 
order V^j's right-hand sides of these equations. The second-order 

approximation is obtained by using the first-order Vco 's in the right-hand 
sides of (3.7-12) and (3.7-13). Finally the values of V_<y[_, and Vyo-p, along 
with the prescribed boundary conditions at each end, are used to evaluate 
(3. 7-14) and (3. 7-15) for k = L and k = T giving the initial and final positions 
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and velocities relative to the launch and target bodies, respectively. All the 
additional equations needed to evaluate the complete solution are given in 
Subsection 3.8. 

The linear solution does not show any explicit dependence on the time 

t = t which appears in the fundamental solution. The solution is a function 

of t however since the constants are evaluated either between t T and t 
o L o 

or between t and t,p (cf. Subsection 3.8). This requires knowledge of not 
only r^(t) but also r^(t) and v_^(t). The time t Q itself is arbitrary to a 
certain extent, and a logical choice would be 

‘o - <‘L + V /2 

The value of r (t ) can then be obtained from the original Lambert solution 

for r (t) and the values of r , (t ) and v , (t ) can be obtained from the funda- 
-o —1' o —1 o 

mental solution with either k = L or k = T and p. = 0. 

If r^(t Q ) is very large then the deviation between the zeroth- and first-order 
solutions may introduce large errors. These errors can be reduced by 
defining a modified linear solution as follows: 

The three position vectors 

L&l) = PL < 4 l> (3.7-16) 

Io (t o } = ^ ( V = £o (t o } + ^^l (t o ) (3.7-17) 

£.Q (t T } = P T (t T ) (3.7-18) 

define the new Lambert problems. The solution of the first Lambert problem 
gives r^(t), shown as the dashed line in Figure 11, and the two zeroth-order 
velocities, v ' (t T ) and v ' (t ). The solution of the second Lambert problem 
gives r^(t), shown as the dotted line in Figure ligand the two zeroth-order 
velocities, y ^(t Q ) and v ^(trp). 
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( 3 . 7 - 20 ) 

( 3 . 7 - 21 ) 

( 3 . 7 - 22 ) 
( 3 . 7 - 23 ) 


where the perturbation terms are obtained from the fundamental solution as 


1! ! I 


?i ^ 


6V^ l - B (t T , t L ) _1 Pf T " + V T '' + p. k T " - A (t T , ^ +Xl + p £' l ) 


+ B» (t T , t o )-vJ(t Q ) - 6 ^ - p 


(3. 7-24) 


6v^ t = c (t T , t T ) + h- i|.) + d (t T , t T ) (6v; T + 6' T + Hi V T ) 


T’ W — L -L 


T’ L/ ' — coL, — L 


- D" (t T , t Q ) v'(t Q ). - 6” -p tIt 1 ! (3.7-25) 

In (3.7-24) and (3. 7-25) the superscripts prime (') and double prime (") 
indicate that the parameter (constant, matrix, function) is evaluated along 
either r^ orr^J respectively. The superscript tilde (~) indicates a special 
partial derivative matrix defined by 

*(t T , t L ) =■*" (t T , t o )®’(t o , t L ) (3.7-26) 

so that 


A(t 


T’ 


V = 


A" (t. 


T’ 


*o> 


A (t 


*l) 


+ B" (t, 


T’ 


t ) C' 
o 


(t 


t L ) 


(3. 2-27) 


B(t T , t L ) = A" (t T , t Q ) B' (t Q , t L ) + B" (t T , t Q ) D r (t Q , t L ) 

(3.7-28) 

C(t T , t L ) = C" (t T , t Q ) A' (t Q , t L ) + D" (t r t Q ) C' (t Q , t L ) 

(3.7-29) 

D(t T , t L ) = C" (t T , t Q ) B' (t o ,t L ) + D" (t T , t Q ) D' (t 0 t L ) (3.7-30) 

The special notation is required because the two zeroth-order solutions, 

r 1 and r", are not continuous in velocity at t = t . Therefore, the transition 
-o -o o 

matrix has discontinuities at t = t . These discontinuities are removed by 

o 

(3.7-26). 
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Equations (3.7-22) through (3.7-30) along with (3. 7-1-4) and (3.7-15) 
constitute a modified linear solution to the boundary value problem. It must 
be evaluated using the same sequence of steps as the linear solution and can 
only be started after r^(t Q ) is found from the linear solution. It thus requires 
more computation time, but forcing the Lambert solutions closer to the 
actual trajectory at t-t Q reduces the size of the perturbation terms with a 
resultant increase in accuracy (cf. the section on numerical results). 

* » - -» r 

The nonlinear interplanetary solution is obtained from the solutions of, a 
sequence of Lambert problems defined by the position vectors 


(t ) 

— o Lt 


El {t L> + p ^- r (t L ) 


(3.7-31). 


(t ) 
-o ' T' 


p T (t T ) + Hi 6.r (t T ) 


(3.7-32) 


where 


6r(t T ) = + I L + [i i 


6r_(t T ) = + V T + p k.' 


(3.7-33) 

(3.7-34) 


Solution of the Lambert problems gives the initial and final velocities, 
V 1 " (t T ) and V’" (t_,). The excess velocities become 

O J -i O 1 


= y » ■ + p y; 


coL 


(3. 7-35) 


-ooT ■" -T + ^ -ooT 


(3.7-36) 


where 


— L = So' (t L> - PL 


y^ = v”’ (t T ) - p T (t T ) 


(3. 7-37) 
(3. 7-38) 
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‘ rgr: 


(3. 7-39) 
(3.7-40) 

This solution also requires a sequence of steps. The zeroth-order solution, 
identicalwith the zeroth-order linear solution, is obtained by putting p = 0 
in (3.7-31), (3.7-32), (3. 7-35), and (3. 7-36). The first-order approximation 
is obtained by putting p = 0 in (3. 7 -33), (3.7-34), (3.7-39), and (3.7-40) 
and using the first order Y^'s in the right-hand sides of (3.7-33) and (3.7-34). 
The second-order approximation is obtained by using the first-order Vco's in 
(3,. 7 -33) and (3.7-34). The nonlinear solutionis shown in Figure 12. 


8 — coL = — L ” ^ — L 

8 -coT = — T “ ^ -T 



Figure 12. Nonlinear Version of Interplanetary Solution 
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It is tempting to define a modified nonlinear solution. by replacing (3.7-33) 
and (3. 7-34) with , . . 

6 r(t L ) = JC’l + + k (3. 7-41) 

6r(t T ) . = + p i” + B" (t T> t Q ) v.', (t Q ) (3. 7-42) ; 


The theoretical basis for making such a substitution is questionable since 


the matrix $ (t^,, t.^), defined by (3.7-26), is not a true transition matrix. 


The nonlinear solution essentially replaces the transition matrix by a new 
Lambert solution to account for nonlinear effects in the partial derivatives. 
The modified nonlinear solution therefore involves the replacement of a - - • 
pseudo-transition matrix with, a Lambert problem -and may not offer any 
improvement over the modified linear solution. " ' • %• - 


All the functions appearing in this section are defined by formulas which 
appear in the next section. 


3.8 FORMULAS FOR THE BOUNDARY VALUE SOLUTIONS 
Each of the preceding solutions requires at least one zeroth-order ellipse 
which is found as the solution of either a standard or a modified Lambert 
problem. Except in the two-impulse moon-to-Earth solution, the zeroth- 
order velocity at t = t^. is used to define a relative velocity 

I k = io (l k ) ■ Ik (t k) (3 - 8 - 1) 


The relative velocity becomes the zeroth-order approximation to the hyper- 
bolic excess velocity, i. e. , 

Y fflk = V k + O (p) (3.8-2) 
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Thus in each of the solutions requiring knowledge of the excess velocity, a 
zeroth-order approximation can be obtained from the Lambert solution. The 
previous sections showed how the zeroth-order approximation becomes the 
starting point for generating successively higher -order approximations to 
V , . The only steps left are to show how and the prescribed boundary 
conditions are used to generate the various parameters appearing in the 
boundary value solutions. The solutions of Subsections 3. 2, 3.3, 3.6, 
and 3. 7 all use the same formulas and these will be given first. The solu- 
tion of :Sub section 3.4 requires a few special formulas to make it fit with the 
other solutions. Finally the solution of Subsection 3. 5 will be considered 
separately. - 

First consider the case where prescribed values of pericenter radius, p^, 
and inclination, i, , must be combined with V_. to determine the inner 
hyperbola. The Cartesian components of are given by the vector 
notation 

2* • <VV *!«> |3 - 8 ' 31 

where the bars indicate inner variables. In addition to the prescribed 
inclination i^, the elements of the hyperbola are (cf. Section B2 of 
Volume 2) 

a k = (3.8-4) 


e k = 1 + P k ^° k 


(3.8-5) 


ctn i. 


cos 


°k 


_ 2 _2 

(u k + V 


v k w k * u k [< u k + v k> ta " 2 s k - w k 


1 / 2 ] 


(3.8-6) 
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ctn i. 


sin ^ = — 


< u k + v k» 


u k w k * v k I (u k + v k> tan 2 ‘k- w k ] 1/2 


(3. 


where 


sin to, 


" ~~[(e^ - 1) 1/2 U,; + Q,. V. 


'k ~k k 


i] 


(3. 


COS to, = 


Sk r (e^. 1) 1/2 V- -Q k U k 


k I u 

e, L 


i] 


(3. 


Uk = U k cos + V k sin ft k 


(3.8 


_ f V k cos S2 k - U k cos 

V k = = ' 
COS 1, 


(3. 8 


w k_ 
sin i. 


(3.8 


and 



(3.8 


8-7) 


8 - 8 ) 


8-9) 


- 10 ) 


-ID 


- 12 ) 


-13) 
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The elements a^., e^, £2^ andu^ are the semimajor axis, eccentricity, 

argument of ascending node, and argument of pericenter, respectively. In 

(3. 8-6) and (3. 8-7) the upper sign is used if the approach to or departure 
th. 

from the k — body is to be over the body, and the lower sign is to be used if 

the approach or departure is under the body. In (3.2;- 12) the upper sign is 

til 

used for departure from the k— body, and the lower sign for approach. 

Additional constants which are used are derived in Sections A12 and B1 of 
Volume 2. They are: 


1 ■ "k = V eok 


(3.8-14) 

A^ = e^ (cos cos £2^ - sin sin £2^ cos i^) 

(3. 8-15) 

( B^ = a^ e^ (cos sin £2^ + sin cos 

£2 k cos i k ) 

(3.8-16) 

Ck “ a k e k s * n W k s * n *k 


(3. 8-17) 

= <\. B k , C k ) 


(3.8-18) 

— ko - — <ok 


(3.8-19) 

— ko = Q k — «k ! n k 


(3. 8-20) 

G ko = Q k log < 2 "k /T k> ^.k^k 

+ ^k 

(3.8-21) 

A k2 = G k — ko /6 


(3. 8-22) 

®kZ =• G k^ko /2 


(3. 8-23) 
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'k2 


°<^ko> ^k2 + G k S ko - 3 ®k2 /2 , 


— k2 


M k(*k Q k*k2< Q kV 

‘ ^k [— k2 " ^k Q k — k2 log ^k + **k Q k — k2] J 


(3.8-25) 


The matrix G(A^ q ) is found from (2. 1-8) with x = A^ q . The matrix is r 
defined as 


G k = G (Bk <*k» 


(3. 2-26) 


The vector (0^/^) I s f° un< I from (2. 2-13) with It must be 

determined by Gaussian quadrature or some other numerical means. An 

* •. ' 
approximate value of A^ can be found from 


— k2 = M k(5k2 lo S 2 l*k - ^ k2 lo « ^k) 


(3.2-27) 


where 


— k2 


G( -ko ) — k2 + H^ko ) -ko — k2 + G k -ko 


/ 2 (3.2-28) 


-k2 


G ^ko> — k2 + «<Ako) G ko ^k2 + G k ^ko • 2 — k2 


(3. 2-29) 


The tensor H( A^ q ) is found from (2. 1-9) with x = A^_ q . 
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Formulas (3.8-4) through (3.8-13) give the constants which completely define 
the inner hyperbola R^ '. Formulas (3.8-14) - (3.8-21) are used to determine 
the behavior of in the overlap domain where is large. The behavior 

of the perturbation, R^, when is large can be found from (3. 8-22) 

through (3.8-29). The key formulas are those for and since they 
appear in (2. 3-5) and eventually in the fundamental solution. 


In the one-impulse moon-to-Earth solution the boundary conditions are not 
given as prescribed pericenter radius and inclination but rather as a pre- 
scribed initial time and position. The initial position Rj^, combined with 
— coM * n (3*4-14), gives the initial velocity Vj^l’ ^he angular momentum 
magnitude, is given hy the relation' ■ ' 


R M1 V coM 


-M 


0 + V*^) sin 


cos 


-Ml 

-acM\ 

R 

V / 

Ml 

COM / 


(3. 8-30) 


■■ - or- 

The inclination is then defined by 


cos 




R u , x 
—Ml 


-Ml\ 


(3. 8-31) 


where e^ is the unit vector in the z direction and x the vector cross product. 
The pericenter radius is obtained from 





1 + V 


coM 





(3. 8-32) 


Thus in the one-impulse solution (3.8-31) and (3.8-32) give the unknown ij^ 
and in terms of and the boundary conditions. These two equations, 

as well as (3.8-30) must be added to the previous set of equations, (3. 8-4) to 
(3.8-29), in order to define the inner hyperbola. 
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The time t^. which can be arbitrarily set ter zero in the solutions of 
Subsections 3.2, 3.3, 3.6, and 3.7, has a fixed value in the one-impulse 
solution. The eccentric anomaly along the inner hyperbola at t - t^ is 
found from 


cosh < F l = >M +R M1> ' (a M; e M> ; 


7(3.8-33) 


s inh F i (Z M1 ’ 5 M i) /( e M a M 1 ) (3.8-34) 


Kepler's equation then gives as 


T M '= (F 1 " e M sinh F l> 7 n M 


(3.8-35) 


i. e. , is the negative of the inner time from pericenter. 


The expressions presented up to this point all pertain to the inner hyperbola 
til 

about the k — body. Another set of formulas defines all of the parameters 
related to the outer solution. Like the inner formulas they involve the 
matrix G defined by (2. 1-8) and the tensor H defined by (2. 1-2), as well as 
the vector function f_ defined by (2-2). These formulas, in order of solution, 
are 




fk = “li'EkV 


N-2 "" 

\S M i [l<£k<V-Pi V ) + i<Ei (t k»] 

i^k 


(3. 8-36) 
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Ek = M k G k Ek “k> 


N-2 


1 M i f G <Pk (t k> - Pi (t k» tec (t k> - Ed (‘k 1 * + ° <Ei <‘ k » Ei <V 

i= 1 - 


i^k 


(3. 


5k = H (Ek (‘k>» 


(3. 


filk = Q k M k-^k> 


(3. 


❖ 


: •• ; • %' .i / 


.^k = -p k Q k G p k /2 


(3. 


[ 5jC ”1 

G k^k-« k -J / 6 - < 3 - 


£4k = Q k M k G ki ( ^k» /6 


<3. 


-10k ( t k' ‘o' " 


/ 


k r 


B(t k' T) ^i < t) + (^TT 


dr 


(3. 


-Ilk (t k’ V 


t 

/ 


D (t k ,T ) Fj (t) - 


—lk 


-2k 


(T-t k )‘ 


- y 


dT 


(3. 


8-37) 

8-38) 

8-39) 

8-40) 

8-41) 

8-42) 

8-43) 
8-44) 
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8-4^ 


8- 4 ^ 




^o 1 

» + * ( v 0 


, t ) + 

) V v ^o 


I*. t ' 

(t v > o 


g-48') 


* &2K 


t 'i r , ^o) 

-\ v c IV ° 


\ T IN 

ft ' V 

S-\^ v ° 

(t ) + & \ ^ 

,, gtv l ° 

+ D vV ° 


IV *o' 


S-SO 1 ) 




g-5^ 


= -3 U* 




xiz) +2W 




- $.\Y 




8-3^ 


* s *V°* 


g<X*> 
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(3.8-54) 


H, 


M k Q k 


i <°k ^ - * & > 


M k Q k 


J k = —2 Ek 


ilk 


k Elk 


— 2k 


°k^ik 


1 _ /~* T * > 

4» n 1 , ( G k c_ 1k - J, a., ) 


3k 


— 4k = G k-lk + »* (G k-lk " J k-lk* 


— 5k ~ °k -lk + H k -lk ‘ ^ J k ~lk 


— 6k = G k &lk + H k ^lk - V- J k ( ^lk + V- ^lk } 


-7k - G k G k -Ik 76 


^8k ‘ 


G k G k ^lk 76 


(3.8-55) 


(3.8-56) 


(3.8-57) 


(3.8-58) 


(3.8-59) 

(3.8-60) 


(3.8-61) 


(3.8-62) 


(3.8-63) 
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-20k (t k* t o ) 


log Q. k (f - t k ) 


K ,„ {t, , t 
— 21k k o 



1 , lo S Q k ( T - t k ) 

4- ib + 4* 7 — 

“ 2 k ( T -V 2 - 3k (T -v 


+ 4* 


4k (t- t k ) 


d t 


(3.8-64) 





D(t k , T )F 2 ( T )- 4lk 


log Q k (t - t k ) 


(T - V 


+ _i_ * log0 k (T - v 

-“(T-t ,) 3 " 3k (r-t ,) 2 


- V 


4k 


(t 


V 


2 ” 5k + 3 — 7k^ 


log Q k (t- t k ) 

' ' < T "~ V 


- +3M 


•6k J -8k 7 (t- t k ) 


dr 


(3.8-65) 
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The formulas (3. 2-36) through (3. 2-52) arise mainly from the expansion of 
the first-order solution, r^, in the overlap domain while those given by 
(3. 2-53) through (3.2-65) arise from a similar expansion of the second-order 
solution, The constants must be evaluated sequentially starting from the 

relative velocity as given by (3.8-1). 


The integral constants K^^., Hi ik* —20k’ anc * —21k mus ^ ke evaluated 
numerically using Gaussian quadrature or some similar technique. This 
requires the integrands to be evaluated at a series of discrete points. The 
first-order integrands are functions of the two-body solution, r^, the partial 
derivative matrices B and D evaluated along r^, and the positions of the N-2 
perturbing bodies as given by an ephemeris, as well as the indicated con- 
stants. The second-order integrands are similar except that they also depend 
on the first-order solution, r^, which is given by (2. 1-13). Since r^ itself 
contains an integral function, the formulas for K^Ok and K ^ ^ are actually 
dotible integrals. This complicates the application of a technique such as 
Gaussian quadrature since each quadrature point of the outer integral must 
be divided into n points for the inner integral. An approximation for r^ which 
reduces and K, ^ to single integrals is discussed in Volume 2. In addi- 

tion, Volume 2 contains explicit formulas for the partial derivative matrices 
A, B, C, and D which must be evaluated to obtain the constants of the outer 
solution. 


The inner constants can be combined to give 


X. 

k 


M, 




, °.k M k e k 

T k 5T k g 2n k 


(3. 8-66) 


while the outer constants can be combined to give 


\ - - ilk <r k - U -iWVV (3.8-67) 

6k = - ilk <‘o -V 1 -5 llk (yy (3.8-68) 
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: Sk * *ik< r k * »' <‘o - *■*&*■%* <‘o - V 1 

• - +t4 k < r k - ‘I -K 20k (‘k- V (3 - 8 - 69 ' 

l k = - ± lk < 2r k 4 l > (t o ' ‘k»' 2/4 - *2k <*o - V /2 

" Atk (t o ' V 4 ±5k ( r k 2 - 2 ) /2 

k±6k« r k +1 l + 3 +7k( 2r k 2 - 5 ) /4 
+ ^3k (3I k ' — 21k (t k* V " ^lk + — k2 

. - M k T k (t lk 4 2 ^ lk - Ek + T M k T k G k I k » « 3 - 8 \ 70 > 

The only other quantities which appear in the boundary value solutions are the 

partial derivative matrices obtained by partitioning the state transition 

matrix, the constants M, defined by (2. 1-7), the dimensionless mass of the 
th ^ 

k — body, p^> and the reference mass, p. As noted previously, many forms 
for the partial derivative matrices are available. Some are in Cartesian 
coordinates and can be used directly. Other forms require coordinate trans- 
formations to obtain the Cartesian expressions. The expressions which are 
given in Volume 2 are in Cartesian coordinates. 

The. two impulse moon-to -Earth solution has not been derived from the funda- 
mental solution and thus requires additional formulas to define the parameters 
appearing in the solution. The hyperbolic solution between the first and sec- 
ond impulses is defined by the following set of orbital elements and angles: 

J= * M (S 1> X V -M (S /) 

I 

N = £ /l £ I 


(3. 8-71) 
(3. 8-72) 


62 


,h = 


V (S + ) 2 


r m< s i> 


a = l/(2h) 


e = (1 + 2h t 2 ) 1/2 


i = co.s * (N • _e^) 


_ N . ei 
sin [2 = — 


sin 1 


cos £2 = — 


N • ^2 


sm l 


+> 


R.. (S VT-V. .(S. 

. , ■=• — M 1 — M 1 

slnhF l = - 1/2- 


cosh Fj = 


a + R M (S 1 } 


a e 


7- a (e 2 -l)^ 2 sinh Fi 
sin £ - — 1 


1 


r m < s i> . 


cos f 


— (e - cos h F i 

"-w- 


[ £3 x N x R m (Sj)] • N. 


sm w = 


R (S 1 ) sin i 


(3. 8-73) 

(3. 8-74) 

(3. 8-75) 
(3. 8-76) 

(3. 8-77) 

(3. 8-78) 

(3. 8-79) 

(3. '8-80) 

(3. 8-8.1) 
(3. 8-82) 

(3. 8-83) 
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£ 3 X -'=-M (S l } 

cos to = — 

R(S ^ ) sin i 



_-3/2 
n = (a) 


f 1 / O 

U = - (Ul sin to + a cos u)/(ae) 

V = (|i?| cos u - a - ^ ^ sin u)/(ae) 

-• • I 

U = U cos £2 - V sin £2 cos i 


(3. 8-84) 


(3. 8-85) 


(3. 8-86) 


(3. 8-87) 


(3. 8-88) 
(3. 8-89) 


(3. 8-90) 
(3. 8-91) 

V = (U,V, W) (3.8-92) 

CO 

The orbital elements, are: a, the semimajor axis; e, the eccentricity; i, the 
inclination; £2, the argument of the ascending node; to, the argument of peri- 
center; and the initial inner time is given by (3. 5-7). These elements and 
the initial time are sufficient to define R. ,-OS) and V..,,-.(S) in a number of 
ways, (2.2-11) being one example. ' 


V = U sin £2 + V cos £2 cos i 


W - V sin i 


The perturbations to hyperbolic motion, which appear in (3. 5-9) and (3.-5- 10) 
require the additional formulas ? 


G.. 

M 

^(Pm^j)) 

(3. 8-93) 

Mm 

= H (Pm^j)) 

(3. 8-94) 

—3 = 

1/24 H m 

(3. 8-95) 
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The outer solution requires the evaluation of four definite integrals given by 


■/ 


—10 (t 2’ t e ) = j B(t 2 , t) FjfTjdr (3.8-96) 


£ll (t 2 ,t e ) 


■/ 


Dft^TtF^TldT 


(3.8-97) 


£ 20 (t 2’ t e ) 


■/ 


B (t 2 , T )F^ 2 ( T )dT 


(3. 8-98) 


r 21 <‘ 2 .v = / D(t 2 > T)F 2 (T)dT 


(3. 8-99) 


These integrals are similar to —20k “““—21k 


and K_ „ which are in the 

— i UK' — 1 IK' — iUK 

fundamental solution. They must be evaluated numerically using a technique 

like Simpson's rule or Gaussian quadrature. The first-order integrands are 

functions of the two-body solution r_ Q obtained from the Lambert solution, the 

partial derivative matrices B and D evaluated along r_ , and the positions of 

o 

the moon and sun. The second-order integrands are similar except that they 
also depend on r^ which, according to (2. 1-13), contains an integral function 
itself. Therefore (3.8-98) and (3.8-99) are actually double integrals and 
require special care in their evaluation. 
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3.9 COMMENTS 

The boundary value solutions presented in the preceding sections satisfy the 
objective stated in the introduction. They are explicit solutions which, 
except for the linear versions of the moon-to -Earth solutions, satisfy the 
boundary conditions exactly without requiring any iterative techniques. Thus 
they offer a distinct advantage over numerical methods which depend on 
iterative techniques to converge to the boundary value solution. It should be 
pointed out however that satisfying the boundary conditions exactly with an 
approximate solution (which.the asymptotic solution does) is not quite the 
same as satisfying the boundary conditions approximately with an exact 
solution (as numerical integration does). Since the asymptotic solution is 
only approximate, an exact solution based on the asymptotic initial conditions 
will most likely not satisfy the terminal boundary conditions exactly. The 

* * * . . . ‘ - ‘ * ' Jt ’’ 

difference between the two methods is the subject of the next part of this 
report. 



Section 4 

" ' NUMERICAL RESULTS 

In order to determine the accuracy of tbe asymptotic boundary value 
solutions, a number of comparisons were made with numerically integrated 
trajectories. Prescribed sets of boundary conditions were used to evaluate 
the asymptotic solutions for Earth-to-moon, moon-to- Earth, and inter- 
planetary applications. Initial conditions determined by the asymptotic solu- 
tions were then used in an N-body numerical integration program. Compari- 
sons of terminal conditions between the asymptotic and numerically 
integrated solutions are used as measures of the accuracy of the asymptotic 
solution. That is, for any terminal condition x the comparisons are given 
as Ax where 

Ax = x (asymptotic) -x (numerical integration) (4-1 ) 

4. 1 EARTH-TO-MOON TRAJECTORIES 

Earth-to-moon trajectories of the type shown in Figures 2, 3, and 4 and 
discussed in Subsections 3. 2 and 3.3 were compared with numerical 
integration. Initial positions relative to the Earth were determined from 
exact solutions of Apollo-type trajectories leaving Earth around February 1, 
1971. 

Five basic trajectories were considered; with variations in the initial 
positions, the total number of trajectories was ten. They cover flight times 
from 60 to 100 hours and inclinations at the moon of -35 and -80 degrees 
(the minus sign indicating an approach under the moon). The boundary 
conditions of all ten examples are given in Table 1. 

(Additional trajectories have also been studied but the major conclusions 
obtained from the results are similar to those obtained from the trajectories 
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BOUNDARY CONDITIONS AND ACCURACIES FOR EARTH- TO- MOON TRAJECTORIES 
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presented here. One of the additional trajectories was a 63. 5 hour transfer 
which was very similar to that used by Carlson (Reference 2). The results 
were in good agreement with his. ) 

The model for the Earth-to-moon trajectories was a 4-body problem (Earth, 
Moon, Sun, Spacecraft) with the positions of the moon and sun determined 
from an analytical ephemeris. 

Since the initial positions for five of the ten trajectories lie close to 
perigee, it was decided to make the comparisons using the nonlinear 
solution, (3. 2-13) through (3. 2-19). The linear solution was used first but 

the results were not satisfactory. The reason for the difference is appar- 

\ 

ently the fact that the linear solution uses the linear partial derivative 
matrices. These are not adequate if one of the endpoints lies in the non- 
linear region of the zeroth-order solution close to perigee. Since the 
nonlinear solution replaces the transition matrix with an exact zeroth-order 
Lambert solution,, the nonlinear effects are included. As the initial position 
moves away from perigee the differences between the two solutions are 
reduced. 

The prescribed initial position _r(t Q ) and the calculated initial velocity 
v(t Q ) were used as initial conditions in the numerical integration program and 
the resulting values of pericynthion time, radius, and inclination compared 
with the prescribed values which are shown on the left side of Table 1. 
Comparisons were made using both the first- and second-order velocities 
obtained from (3.2-16) and are shown on the right side of Table 1. The 
results can be divided into groups which show various trends in the accuracy 
of the asymptotic solution. The data for each group are presented in the 
following format. 

Group Case Order of 

Designation Number Solution 


At p M Ap M Al M 
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The groups are: 

A. 101 

1 

11 

-18 

-0. 05 

102 

1 

36 

-92 

-0. 05 

103 

1 

45 

-210 

-0. 11 

In this group of three trajectories, 

the variable boundary 

condition is 

the 

time of flight, and the results show that first- 

order accuracy is degraded as 

the time of flight increases. 
B. 101 

2 

2 

-6 

- 0. 15. 

102 

2 

-8 

-299 

t 0 , 3 3 

103 

o 

2 

-154. . 

-4497 

-2. 12; 


In this group, the variable boundary condition is again time of flight but the^. - 
results are now second order. The same degradation of accuracy as in A is 
apparent but the degree of degradation is much more marked. For No. 101 
the second order is better than the first, but for No. 103, the first order is , 
better. 

These results indicate that whatever is degrading the first-order accuracy 
is having a highly adverse effect on the second-order solution. The cause 
is most likely the deviation of the first-order solution from the zeroth-order , 
solution. When this deviation is large, as it is in the long -flight-time 
trajectories, it evidently causes the asymptotic solution to diverge before 
reaching the second-order term. Thus the second-order error is larger 
than the first. 


C. 

101 

1 

11 

-1.8 

-0.05 


104 

1 

1 1 

-21 

-0. 04 


101 

2 

2 

-6 

0. 15 


104 

2 

3 . 

4 

0. 08 


This group shows the effect of choosing the initial position to be a point one 
hour out on trajectory No. 101, thus giving No. 104. By delaying the initial 
position, it occurs at a true anomaly of 118 rather than 6 degrees and at a 
radius of almost 13, 000. rather than 3, 544 nautical miles. This delay has 
little effect on either first- or second-order accuracy. 
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103 

1 

45 

-210 

0. 11 

106 

1 

47 

-162 

-0. 34 

103 

2 

-154 

-4497 

-2. 12 

105 

2 

-12 

-699 

-0. 65 

106 

2 

1 

-421 

-0. 45 


This group is similar to C except that the time of flight is near 100 hours 
rather than 60. A delay of one hour (No. 105) along No. 103 increases the 
true anomaly by 106 degrees and the radius by 9, 360 nautical miles. A delay 
of two hours (No. 106) increases the true anomaly another 15 degrees and 
radius another 8, 264 nautical miles. The two -hour delay has very little 
effect on the first order results but the one and two hour delays have a 
marked effect on the second order results reducing the time of flight error 
from 154 minutes down to one minute and the radius error from almost 
4, 500 nmi down to 421 nmi. 


Carlson (Reference 2) has showii such a trend in the first-order linear 
results for a 63, 5-hour trajectory (starting at perigee). He shows steps 
(delays) of 1, 3. 5, 8, 5, 13. 5, 23. 5, 33. 5, . . . hours. Accuracy is improved 
with each step out to 8. 5 hours after which it slowly begins to degrade. He 
found no such trend using a first-order nonlinear solution, which agrees with 
the first-order comparisons shown here and in C. 

The improvement in accuracy shown here as the initial position is moved 
away from the Earth is evidently due to a decrease in the deviation between 
the zeroth- and first-order solutions. This deviation is small for No. 101 
thus a delay of one hour has little effect. For No. 103 however, the deviation 
is initially quite large and delays of one and two hours decrease the deviation 
and increase the accuracy. The deviation of the solutions can be measured 
by the magnitude of the first order 6v(t Q ). The magnitudes are 


No. 

101 

104 
103 

105 

106 


Delay (hr) 
0 
1 
0 
1 
2 


6v(t 0 ) (ft/sec) 

149 

54 

845 

254 

210 
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Since neither 101 nor 104 have excessively large values of 6v(to) the second 
order solutions are not adversely affected. The large value for No. 103, 
however, results in large second order errors. 

E. 107 1 11 2 -0.65 

108 1 46 -94 -4.24 

This group is similar to A except the inclination at the moon is -80 rather 

than -35 degrees and data for the 80-hour trajectory has not been included. 

The same degradation with increasing time of flight is apparent and, except 
for inclination, the errors are similar in magnitude to those of A. 

F. 107 2 -64 -644 -8.86 

108 2 -113 -2499 -25.90 

This group is similar to B except for the difference in inclination and the 

lack of an 80 hour trajectory. It again shows the degradation of the second- 
order results as flight time increases but, except for inclination, the degrada- 
tion is not as marked as in B. 


Another trend is also apparent when B and F are considered together as 
follows: 


G. 

101 

2 

2 

-6 

0. 15 


107 

2 

-64 

-644 

-8. 86 


, 103 

2 

-154 

-4497 

-2. 12 


108 

2 

-113 

-2499 

-25. 90 


Comparison of Cases 101 and 107 shows that the second-order accuracy is 
degraded as the inclination goes from -35 to -80 degrees. Comparison of 
cases 103 and 108 shows that, except for inclination, the second-order 
accuracy is improved as the inclination changes. (The anomalous behavior 
of the inclination cannot be explained. ) 

Thus inclination at the moon has an effect on the second -order accuracy, the 
nature of the effect being dependent on the time of flight. No such trend is 
apparent in the first-order results. 
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H. 

107 

1 

11 

2 . 

-0. 65 


109 or 

1 

11 

r -11- 

-0. 48 


107 . 

2 

-64 

-644 

-8. 86 


109 

2 

1 

-16 

0. 36 


This group is similar to C except for the difference in inclination. No. 109 
is obtained by delaying the starting position one hour along 107. Again the 
delay has little effect on the first-order results but results in a marked 
improvement in the second-order results where the errors were originally 


quite large. 

The reason 

for the improvement 

is discussed in D. 

- 

I. 

108 

1 

46 -94 

-4. 24 

1 

110 

1 

47 -89 

-3. 07 


.108 

2 -113 -2949 

-25. 90 


no 

2 

3 -353 

-7. 71 


This group is similar to H except that the time of flight is near 100 hours 
rather than 60. The trend is identical to H. A delay of two hours has little 
effect on the first-order results but a significant effect on the second-order 
results. The reason is again the same: an initial position further from the 

Earth reduces the deviation between the zeroth- and first-order solutions and 
thus increases the accuracy. 

Several observations can be drawn from A through I and from Table 1 as a 
whole. The first would be that the second-order solution is much more 
sensitive to boundary conditions than is the first order. For instance first- 
order pericynthion radius errors range from -210 to +2 nmi while the 
second-order errors range from -4, 497 to +4 nmi. In some cases the 
second order is better than the first; in other cases the first order is better. 
The point of minimum error for the asymptotic solution is evidently less 
than second order in some cases, and adding the second-order terms causes 
the solution to diverge. {In general, an asymptotic expansion in powers of p 
will most closely approximate the function which it represents after n terms, 
and then diverges as n is increased. The optimum value of n is a function of 
p and other parameters such as boundary conditions. It is difficult to 
determine a priori. ) 
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The second observation is that second-ord'er errors are increased by 
starting the trajectory close to perigee, by increasing the time of flight, and 
by choosing certain values of inclination at the moon. Each of these effects 
can be related to the size of the first-order velocity correction 


6v(to) = |V (ty - lo^oJ ^ 


(4.1-2) 


When the magnitude 6v(t 0 )'is large the error is large and when the magnitude 
is small the error is small. When v(t Q j is large, the first-order solution 

_r , tends to grow rapidly with time; From (2.1-6) it can be seen that the 

1 2 

force is' a function of r^ and this results in large values of the second- 
order integrals K_ 2 q M and ^ ^ Since 


- r l - (£ - £o)/hl ' (4.1. 3) 

its effect could be reduced by choosing to be closer to the actual 
trajectory _r. - However, the condition dictated by (2. 2 -4) prohibits from 
being chosen arbitrarily, i. e. , it must pass through the center of the moon. 

If (2. 2-4) is relaxed then it might be possible to choose an r_ Q which would 
minimize (4. 1-3) but such a step causes rather severe complications in the 
analytical derivation of the solution and is not well enough understood at 
this point to make any modifications along this line. (For a further discus- 
sion -of (2; 2-4) see Section A10 and the first parts of Sections All and A17 
in Volume 2. j 

• ■ . ■» ■ , " 

The magnitude of 6v(t 0 ) should be proportional to the perturbations 

experienced over the entire length of the trajectory. Numerically integrated 

_5 

trajectories show a variation in the orbital elements of less than 10 over 
the first two hours of flight after leaving perigee, indicating that the moon 
and sun perturbations are small over this interval. The large changes in 
the magnitude of, 6v(t 0 ) as the starting point moves away from perigee (cf. , D)' 
must therefore be attributed to something other than the effect of the moon 
and the sun. The evidence indicates that the cause is the large difference 
between r and r + u. r., i. e. , a fictitious force arises when r is a poor 
approximation to the actual solution. 


74 


Since a delay in the starting position to one or two hours after perigee never 
degrades the accuracy bu;t,.tends to improve it when the v error.s are large, 
and since the effect of the real perturbing forces is small over the initial 
one- or two-hour interval it may be possible to construct a more uniformly: 
valid solution by applying either a two-body or a perturbed two-body solution 
out to a true anomaly of 120 to 135 degrees and then applying the asymptotic 
boundary value solution from there to the moon. Certain modifications would 
be required to make the velocity continuous between the two solutions, but 
such changes are relatively minor. Lack of time prohibited making such a 
modification in this study although a somewhat related solution was 
incorporated into the interplanetary solution. It is termed the modified 
linear solution in Subsection 3. 7 and the results are shown in Subsection 4.-4. 

A. third observation is that since accuracy is increased as the starting point 
is moved away from perigee, the midcourse application is probably better 
than the full Earth-to-moon application. This is especially true if a simple 
form of the asymptotic solution is desired, since it has been amply pointed 
out that starting from perigee probably will require a composite solution to 
obtain uniform accuracy. 

The fourth and final observation is that starting at perigee will always 
involve rather large sensitivities to initial errors. The initial velocity . 
calculated, even without any modification, may differ by only a small amount 
from the actual velocity which satisfies the boundary value problem exactly. 
Whether or not these small differences are tolerable depends on the 
particular use to which the solution is being applied, 

4. 2 TWO-IMPULSE MOON- TO - EARTH TRAJECTORIES 

Two -impulse moon-to- Earth trajectories of the type shown in Figures .7 and 
8 and discussed in Subsection 3. 5 were also compared with numerical 
integration. The initial position and velocity were taken from a patched 
conic optimum 4-impulse solution obtained from NASA-MSC. The initial 
time was zero hours on February 8, 1980 (Julian Date 2444277. 5) and the 
initial position and velocity were _ . 

—Ml = (263. 24, -923.75, 272, 72) nautical miles (4,2-1) 

V' = (-1125.44, 1181.10, 5086. 90) feet/second (4. 2-2) 
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This corresponds to a 60 -nmi circular orbit about the moon with an inclina- 
tion of 97. 582 degrees. The first impulse was chosen to match the patched 
conic solution by putting 

l l = 0.45805 (4.2-3) 

This resulted in a hyperbolic orbit after the first impulse with a pericenter 
radius of 998, 49 nmi, an eccentricity of 1. 1259, and the same inclination 
as before. 

Equations (3. 5-8) and (3, 5-9) were used to determine the position and 
velocity at five-hour intervals out to 25 hours. The magnitudes of the two- 
body position and velocity, R Mq and V Mo , the position and velocity perturba- 
tion, 6R^, and 6V, ,, and the errors when compared to numerical integration, 
M M 

| 6 R M | and I^V^I . are shown in Table 2. The results show the solution 
defined by (3. 5-8) and (3. 5-9) to be quite accurate with position and velocity 
errors running 0. 1 percent or less. It is interesting to note that at 20 hours 
the predicted radius was approximately the radius of the moon's sphere of 
influence and the perturbations 6R M and &V M indicate the error of a patched 
conic solution at this point. 


Table 2 

POSITION AND VELOCITY MAGNITUDES AND 
ERRORS OF PERTURBED HYPERBOLA 


t 2" t l 

(hr) 

r mo 

(nmi) 

6R, , 
— M 

(nmi) 

im m i 

(nmi) 

V M0 

(fps) 

6V M 

(fps) 

I- v mI 

(fps) 

5 

10, 842 

5. 7 

0. 9 

2,975 

4. 9 

0.4 

10 

18, 797 

39. 6 

5. 1 

2, 574 

17.4 

0. 8 

15 

26, 070 

124.0 

13. 2 

2,404 

37. 2 

1. 1 

20 

32,992 

281.4 

24. 9 

2, 307 

64.6 

1. 0 

25 

39,694 

535. 2 

39. 7 

2,243 

100. 0 

2. 5 
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At each of the five points in Table 2, a second impulse was calculated using 
the linear solution (3. 5-19). Entry times were chosen to make the total 
flight time, t - t^, equal to 120, 100, and 80 hours. In addition, an 
impulse was added at t^ - t^ = 5 hours which resulted in a 60 -hour flight 

time. (Flights of 60 hours flight time from t^ - t^ = 10, 15, 20, and 25 
hours were not possible without going to hyperbolic transfers. This was. 
because t - t~ for each of these points was less than the parabolic flight 
time from that point to Earth. ) These trajectories are summarized in 
Table 3, They are grouped according to the time of the second impulse in 
order of decreasing total flight time. Magnitudes of the second velocity 
impulse are shown, the magnitude of the first impulse was the same for 
all cases, i. e. , 

AV(t.) = 2, 447 fps (4.2-4) 

The Earth entry conditions were 

r = 3, 594 nmi (4.2-5) 

e 

Ye = 0. 0 deg (4.2-6) 

i = 30. 0 deg (4.2-7) 

The choice of zero for Y g meant that entry coincided with perigee. Errors in 
perigee time, radius, and inclination are shown as well as errors in entry 
time and flight path angle for those trajectories where the numerically 
integrated solution passed through the entry radius, i, e. , had positive perigee 
radius errors meaning the numerically integrated perigee was below the 
perigee (entry) radius. 

The numerical integration program was set to cut off at flight times of 125 
hours. In the cases where this limit was reached, the radius at that time 
is shown to give an indication of how far off the trajectory was. 
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TIME INTERVALS, VELOCITY IMPULSES, AND ACCURACIES FOR 
TWO- IMPULSE MOON- TO- EARTH TRAJECTORIES (Page 1 of 2) 
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’’-'Numerically inlegrated perigee radius is greater than prescribed entry radius and no entry occurs. 



Table 3 

time intervals, velocity impulses, and accuracies for 

TWO- IMPULSE MOON- TO- EARTH TRAJECTORIES ( Page 2 of 2) 
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❖^Patched conic solution obtained from NASA-MSC. 



The results shown are second-order results except for the 120-hour flights 
and the flights where the second impulse occurred at the sphere of influence. 
In these cases zeroth, first- and second-order results are all shown. (The 
order of the solution applies only to the phase after the second impulse. 

The hyperbolic phase used the perturbed solution in all cases since the 
accuracy of this phase had already been determined. ) 

The results show three basic trends. First, as t^ - t^ increases with .. 

t g - t^ held fixed, the accuracy improves. This is to be expected. As the 
point of the second impulse moves away from the moon, the outer solution 
which is used to derive (3. 5-19) becomes a better approximation to the actual 
solution. If the initial point is too close to the moon (a distance of order 
then the outer solution must be replaced by a matched solution obtained from 
the fundamental solution. Since the primary purpose of the two-impulse 
solution is to study trajectories where the second impulse occurs away from 
the moon, the derivation was made accordingly. Thus values of t^ - tj of 5 
and even 10 hours represent marginal applications for this type of solution, 
and the errors are not entirely unexpected. 

The second trend occurs when t 0 - t, is held fixed and t - t, decreased. 

The trend in each group is toward improved accuracy. This should also be 
expected if the Earth-to-moon results of the previous section are considered. 
It was shown there that increasing flight time degraded accuracy. The same 
trend appears here with decreasing flight time giving increased accuracy. 

The cause is probably the same as in the Earth-to-moon cases, i. e. , the 
deviation of the first-order solution from the zeroth-order solution increases 
with increasing flight time and adversely affects the accuracy. A second 
cause may be the use of the linear rather than the nonlinear solution. - All of 
the post-second impulse trajectories had transfer angles between 170 and 
180 degrees, the 120-hour trajectories all being around 178 degrees. The 
linear partial derivative matrices may introduce considerable error in such 
cases since both end points lie in highly nonlinear regions of the zeroth-order 
solution. Since the nonlinear solution does not use the partial derivative 
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matrices (except in the definite integrals), it could be expected to give better 
results. Time did not permit an evaluation of the nonlinear solution. 

The third trend apparent is that the first-order solution gives better results 
than the second order in all but one of those cases where the zeroth-, first-, 
and second-order solutions have all been compared with numerical integra- 
tion. However, a comparison of Cases 201, 205, 208, 211, and 214 shows 
that the difference between the first- and second-order solutions is decreas- 
ing as t^ - t^ increases and finally when - t^ equals 25 hours the second 
order is better. This also can be attributed to the size of the deviation 
between the zeroth- and first-order solutions. When this deviation is large, 
as it is when the second impulse occurs close to the moon, then the second- 
order accuracy is adversely affected. As the deviation decreases, the 
adverse effect is diminished and the second-order accuracy is improved. 

The cases where t^ - t^ equals 20 hours, Cases 211 to 213, are interesting 
since the zeroth-order results are nearly equal to what a patched conic 
solution would give. Prior to the second impulse, the solution is a perturbed 
hyperbola, but the perturbations are not excessively large and the solution 
is close to a pure two-body solution. After the second impulse, the zeroth- 
order solution is a two-body ellipse. Since the second impulse occurs at 
the sphere of influence, the zeroth-order solution is nearly patched conic 
(actually slightly more accurate). The results show that both the first- and 
second-order solutions are an order of magnitude more accurate than the 
zeroth-order solution and that the zeroth-order solution predicts a velocity 
impulse which is less than that actually required to satisfy the boundary 
conditions. This would indicate that patched conic solutions may give 
considerable error if used to determine the optimum time for adding the 
second impulse. 

This point is emphasized by considering No. 217. The errors for this case 
were obtained by numerically integrating from the initial conditions predicted 
by the patched conic program used at NASA-MSC and adding the indicated 
impulse at 6. 3 hours. Although the solution is supposedly optimum, the 
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terminal errors are quite large and the second impulse would probably have 
to be increased by several hundred feet per second to cancel out the errors. 
The magnitude of the second impulse would then approach that of No. 209 
which was the minimum of all the second-order impulses calculated by the 
asymptotic solution. However, for No. 209 the second impulse occurs 
8. 7 hours later and the total flight time is 20 hours less than No. 217, 
indicating quite a change in the possible optimum conditions. 

A final observation is that the errors in time of flight, perigee radius and 
entry flight path angle shown in Table 3 are never as small as one might 
expect from a second-order theory. As discussed previously, the overall 
accuracy could be expected to improve by going to a nonlinear solution, but 
problems pertaining to sensitivity might still remain. Thus, as in the Earth- 
to-moon solution, the intended use of the solution will determine to a large 
extent what accuracies are acceptable. If the accuracies shown here, 
particularly those for long flight times, are unacceptable then further 
numerical analysis would be warranted to see if the errors could be reduced 
by either using the nonlinear solution or by formulating a composite solution 
in which the deviations from the zeroth-order solution are less. This latter 
approach would be similar to the combined two -body /asymptotic solution 
suggested for Earth-to-moon trajectories. 

4. 3 INTERPLANETARY MIDCOURSE TRAJECTORIES 

Interplanetary trajectories of the type shown in Figure 9 and discussed in 
Subsection 3. 6 were compared with numerical integration. Two Earth-to- 
Mars reference trajectories were chosen; the first was a 244-day transfer 
leaving Earth on November 11, 1964 and the second a 184-day transfer 
leaving Earth on May 19, 1971. The first reference trajectory was used by 
Carlson (Reference 2) and is presented here for comparison with his 
results. The second reference trajectory is similar to that actually flown 
by the 1971 Mariner mission. 

The prescribed boundary conditions at Mars are shown in Table 4. The 
initial positions were determined along two-body solutions intersecting 
Earth and Mars at the departure and arrival dates. The initial velocity was 
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Table 4 


TERMINAL BOUNDARY CONDITIONS FOR 
MIDCOURSE-TO-MARS TRAJECTORIES 


Reference 

Trajectory 

t 

PT 

(Date) 

(Hour) 

PT 

(nmi) 

*T 

(deg) 

1964—244 Day 

July 13, 1965 

11:59:45 

1, 898 

33. 9 

1971-184 Day 

Nov 19, 1971 

12:00 

2, 000 

30. 0 


then calculated using the linear solution. The position and velocity were 
used as initial conditions in the numerical integration program and the 
resulting boundary conditions compared with the predicted values. The 
comparisons are shown in Table 5, Cases 301 to 306 corresponding to initial 
positions along the 1964 244-day trajectory, and Cases 307 and 308 corre- 
sponding to initial positions along the 1971 184-day trajectory. The time to 
go to pericenter and the distance from Mars are shown for each initial posi- 
tion as well as the differences (errors) in pericenter time, radius, and 
inclination. 

The model used was a seven-body problem with the positions of Venus, 
Earth, Mars, Jupiter, and Saturn obtained from an analytical ephemeris. 

The 1964 reference trajectory is a relatively low-energy Earth-to-Mars 
transfer with a heliocentric transfer angle of about 178 degrees. The six 
midcourse points, 40 days apart, correspond to some of the initial positions 
investigated by Carlson (Reference 2). In each case both first- and second- 
order results are shown. 

The first-order comparisons agree very well with Carlson's first-order 
results. They show an increasing degree of accuracy from 220 to 100 days 
and then a trend toward decreasing accuracy as the time of flight is further 
reduced below 100 days. The large error at 220 days can be attributed to 
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Table 5 

ACCURACIES FOR MIDCOURSE- TO- MARS TRAJECTORIES 
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several causes but the most probable is that the starting point is in a region 
relatively close to the Earth where the midcourse solution may not be 
applicable. That is, the assumptions used to derive the midcourse solution 
are not satisfied if one of the perturbing bodies other than the target body is 
relatively close to the starting point. If the starting point is within a distance 
of order p of such a perturbing body then the interplanetary solution of 
Section 3.7 must be used. The starting points for Cases 301 and 302 
'evidently lie between the domain of validity of the two solutions. The 
increasing error as tp’p-to goes below 100 days can be attributed to certain 
(terms proportional to (tprp-t 0 )“* which are ignored in the first-order solution. 


The second-order solution shows a definite improvement in accuracy for 
Cases 301 and 306. InCases 302 to 305 the second-order shows an improve- 
ment in time of pericenter passage but the results for radius and inclination 
are inconsistent, sometimes better and sometimes worse than the first 
order. The mixed results for radius and inclination are probably due to the 
fact that both the first- and second-order values may be within the accuracy 
jlimits of the numerical integration itself, 
j 

The 1971 reference trajectory is a somewhat higher-energy Earth-to-Mars 
transfer with a heliocentric transfer angle of about 142 degrees. The two 
midcourse points occur two and four months after launch and lie in the 
region where the midcourse solution can be expected to work well. For each 
case zeroth- , first- , and second-order results are shown. 


The zeroth order solution, which is identical to the massless planet conic 
approximation often used for interplanetary trajectory analysis, results in 
'rather large errors. These errors are reduced significantly by the first 
j and second order solutions and, except for the pericenter radius error in 
- Case 307, the second order is better than the first. 

The magnitude of the first- and second-order errors in Cases 307 and 308 is 
difficult to explain. Since the 1971 trajectory is a higher energy one (i. e. , has 
1 less flight time) than the 1964 trajectory, the results from the lunar trajecto- 
ries would indicate that the errors should be less. Also the two points are at 
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distances from Mars which are comparable to Cases 304 to 306. Yet both 
the first- and second-order errors in Cases 307 and 308 are considerably 
larger than those in Cases 304 to 306. These differences are evidently due 
to the combination of time of flight, pericenter radius, and inclination pres- 
cribed for each of the reference trajectories, but a detailed inspection of each 
of the cases does not give any hint as to why this should be so. 

The final observation regarding Cases 307 and 308 Ls that the first- and 
second-order results show improvement as the time to go decreases. This 
is in agreement with Cases 301 to 306. 

The interplanetary midcourse solution appears to be an excellent application 
for the asymptotic solution. There are two basic reasons for this. The first 
is the magnitude of the small parameter p. In the lunar cases p is 10 
while in the midcourse-to-Mars application pis 10“^. This causes the 
asymptotic expansions to converge much more rapidly. 

The second reason is less apparent from a theoretical point of view but just 
as important as the size of p: it is the fact that the location of the initial 
position does not have as strong an effect on accuracy as in the previous 
results. In the lunar examples the zeroth-order solution has an eccentricity 
close to unity and the region close to perigee (and apogee) has a highly non- 
linear behavior. In interplanetary applications the corresponding eccentri- 
cities are 0.25 or less and the initial position, even if close to perihelion, 
does not lie in a highly nonlinear region. 

I , f 

The combination of small p and nearly circular zeroth-order outer solutions 
causes the deviation between zeroth- and first-order solutions to be rela- 
tively small when both start from a common position as in the midcourse 
solution. The results of the next section show that this is no longer true 
when the initial position is close to another body. 
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4.4 INTERPLANETARY TRAJECTORIES 

Interplanetary trajectories of the type shown in Figure 10 and discussed in 
Subsection 3. 7 were also compared with numerical integration. The two 
reference trajectories are the 1964 and 1971 examples of the previous 
section except here they are considered in reverse order. The same seven- 
body model was used. 

The prescribed boundary conditions are shown in Table 6. The conditions 
at Earth are typical of departure from a low-Earth orbit. In Cases 401 and 
402 the terminal conditions are typical of a close approach, while in Cases 403 
and 404 the pericenter radius at Mars is rather large. This value was 
chosen since a numerically integrated solution with identical boundary condi- 
tions was available for comparison. 

' 1 

For each case, the initial position and velocity at the Earth were calculated 
from (3. 7-14) and (3. 7-15) using either the linear or modified linear solution. 
The position and velocity were then numerically integrated up to a close 
approach at Mars and the resulting boundary conditions compared with the 
prescribed values. The comparisons are shown in Table 7. 

Case 401 shows zeroth- , first- , and second-order results of applying the 
linear solution to the 1971 trajectory. The results show the expected 
improvement going from zeroth to first order but then a significant degrada- 
tion in accuracy going from first to second order. 

i 

Case 402 shows the results of applying the modified linear solution. The 

_3 

midpoint was offset by almost 75, 000 nautical miles or 10 in dimensionless 
units. The results show that except for zeroth-order inclination and second- 
order time of pericenter passage, the modified, solution caused a slight 
degradation in accuracy. 

Since the zeroth- and second-order solutions showed such large errors 
for the 1971 trajectory, only the first-order solution was compared with 
numerical integration for the 1964 trajectory. The first-order results of 


o 


87 



INITIAL AND TERMINAL CONDITIONS FOR EARTH- TO- MARS TRAJECTORIES 
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Table 7 

ACCURACIES FOR EARTH- TO- MARS TRAJECTORIES 


Trajectory 

Number 

Heliocentric 
Transfer 
Angle (deg) 

Midpoint 

Shift 

(nmi) 

Order 

At 

PT 

(hr) 

Ap t 

(nmi) 

Ar T 

(deg) 

401 

142.4 


0 

101. 4 

-71, 702 

65. 4 




1 

8. 0 

-8, 289 

-0. 5 




2 

472. 5 

-32, 370 

-9. 9 

402 

142.4 

74, 642 

0 

1.17. 5 

-83, 357 

13. 5 




1 

9. 0 

-8, 819 

-0. 8 




2 

456. 7 

-32, 626 

-9. 9 

403 

177. 6 


1 

8. .7 

-22, 545 

0. 8 

404 

177. 6 

920, 987 

1 

3. 6 

-9, 005 

-2. 4 


both the linear and modified solutions are given by Cases 403 and 404, The 
modified solution resulted in a reduction of over 50 percent in the time and 
radius errors but caused a slight increase in the already small inclination 
error. The reason that the modified solution gave better results in this 
example is probably due to the size of the midpoint shift which is 1 0" c in 
dimensionless units. Since the midpoint shift was an order of magnitude 
larger than in the previous example it could be concluded that the linear 
solution (403) contained significantly more error (than 401) due to the larger 
deviation from the zeroth-order solution. 

It should be pointed out that a solution of any given order has a minimum 
error associated with that order. Solutions like the modified linear or the 
nonlinear are only attempts to reduce the error to its minimum value. If 
the linear solution itself is close to the minimum error then the other solu- 
tions will not result in any improvement. This is evidently what happened 
in Case 402. 
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The second-order errors are much larger than anticipated and are not removed 
by going to the modified solution. It is felt that the cause of the large error 
is twofold. First, as shown in Figure 10, the zeroth-order solutions do not 
go through the true midpoint in either the linear or modified linear solutions. 
Therefore there is always a finite difference between zeroth and first order 
and it has been shown in previous sections that this can cause rather large 
second-order errors if the deviation becomes large. 

In the lunar and interplanetary midcourse solutions, the first-order perturba- 
tion, r^, always vanishes at t = to and this helps to reduce the maximum 
value which rj may attain at some later time. In the interplanetary solution 
such is not the case and the large values which rj may attain cause 
excessively large deviations from zeroth order. In order to illustrate this 
point consider the difference between the zeroth- and first-order midpoints 
for Cases 401 and 403. These differences are equal to the midpoint shifts 
in Cases 402 and 404, i.e., 10~3 and 10"^. But according to (3.7-17) the 
shift in the midpoint should be order por 10“^. (If the earth is the reference 
body then jj would be 10"^. In either case the midpoint shift should be the 
same, theoretically between 10“^ and 10“^.) Even if these differences are 
reduced by one or two orders of magnitude they still remain larger than 
order p and therefore introduce error. 

The second cause probably lies in the evaluation of the integral K 21 L which 
eventually ends up in q In a H other solutions which use this integral it is 
not used in calculating the initial position or velocity. That is, the value of 
the initial velocity perturbation 6v(to) is a function of the first order value of 
Vaa^. and the constant Tk use & only in obtaining the second order value of 
Voo^. This means that q^- is used only in calculating non-prescribed boundary 
conditions at the terminal end of the trajectory. In the interplanetary solution 
solution, however, q ^ is used to determine the second order value of 
[cf. (3. 7-10) and (3. 7-12)] which in turn is used to determine second order 
orbital elements and second order position and velocity at the launch planet 
(cf. (3. 7-14), (3. 7- 1 5) and Section 3. 8). Therefore the importance of the 
integral increases significantly in the interplanetary solution. 
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Of all the definite integrals which must be evaluated IS21k most 

difficult because it has the most complex integrand. Since the integrals 

cannot be evaluated exactly but only by some numerical approximation they 

always include some error. If they are large then the error is correspond- 

8 

ingly large. In Case 402 the magnitude of ^ was 10 . A one percent 
error in the Gaussian quadrature would be 10^. The effect of this error can 
be determined by multiplying by since I$21I_, is a second-order integral. 

The net effect is 10"® or almost order p| If attempts are made to reduce 
the error by enough orders of magnitude so that it is no longer a problem, 
then the computation time increases considerably. 

The large second-order errors indicate that a better zeroth-order solution is 
needed for this application. The first-order solution however, works well 
and offers considerable improvement over the massless-planet (zeroth-order) 
approximation. Even though the first-order errors in Cases 401 to 404 may 
appear somewhat large, the actual initial positions and velocities are quite 
close to those needed to satisfy the boundary conditions with an exact 
solution. The errors between the asymptotic solution and numerical 

O' ... 

integration when both satisfy the boundary conditions are shown for 
Cases 403 and 404 in Table 8. Notice that in the linear solution (403) the 
first-order errors are large and the second order actually slightly better. 

By going to the modified solution (404) the zeroth- and first-order errors 
are reduced significantly while the second-order errors are reduced only 
slightly. 

The first-order error for Case 404 dramatically illustrates the problem of 
sensitivity which has been mentioned earlier. Errors of 14 nmi in position 
and 77 fps in velocity (out of a total velocity of 37, 858 fps) result in errors 
of 3.6 hours and 9, 005 nautical miles at closest approach to Mars. Since 

i. ' 

problems of sensitivity cannot be eliminated, the errors shown for first 
order in Table 8 are probably acceptable. 

4. 5 COMPUTATION TIMES 

The asymptotic solutions require three basic types of calculations: (1) Lam- 
bert solutions, (2) Gaussian quadrature, and (3) general calculations 
including the boundary value solution equations, partial derivative matrices, 
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Table 8 


INITIAL POSITION AND VELOCITY ERRORS 
FOR EARTH- TO- MARS TRAJECTORIES 


Trajectory 

Number 

Order 

(nmi) 

l— piJ 

(fps) 

403 

0 

601 

2, 706 


1 

84 

460 


2 

47 

321 

404 

0 

96 

561 


1 

14 

77 


2 

35 

326 


conversions to dimensionless numbers, etc. The first two require numerical 
techniques (including iterations in the Lambert solutions); the third involves 
straightforward evaluation of explicit expressions. The calculations can be 
accomplished quite rapidly on a high-speed computer. 

FORTRAN programs which evaluate the asymptotic solutions were run on a 
CDC 6500 computer as was the numerical integration program. In the 
asymptotic solutions the Lambert problems required' approximately 0. 3 
seconds per solution while a 60-point Gaussian quadrature routine required 
approximately 0.4 seconds for each of the N-2 perturbing bodies. The total 
computation times for each of the different solutions as well as the corres- 
ponding numerical integration times are shown in Table 9. 

In the Earth-to-moon solution a large part of the time is spent in Lambert 
solutions, while in the interplanetary solutions a considerable portion of the 
time is spent calculating the effect of the perturbations through the Gaussian 
quadrature routine. Since the interplanetary solution is derived from two 
fundamental solutions rather than one, it requires almost twice as much time 
as the midcourse solution. The modified linear solution adds two additional 
Lambert solutions plus an extra first-order solution which pushes the time 
over 10 seconds. 
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Table 9 

COMPUTATION TIMES ON CDC 6500 COMPUTER 
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Whereas, the, asymptotic solutions are boundary value solutions, the numerical 
integration, runs with which they were compared are only initial value 
solutions based on the. asymptotic initial conditions. Even so, the computa- 
tion times for the asymptotic solutions are from 6 to 38 times faster than 
those for numerical integration. In order to obtain a boundary value solution 
from numerical integration, a. hunting or, search procedure must be added 
which significantly increases the .computation time. If the asymptotic 
solution is then compared with the numerically integrated boundary value 
solutions, the computation times are 25 to. 150 times faster. An actual com- 
parison is shown in Table 9 for interplanetary midcourse trajectories. 

It should be noted that computation time is significantly affected by the 
number of quadrature points needed to evaluate the constants. A 120-day 
interplanetary midcourse Solution .was evaluated with 12, 36, 60, and 96 
quadrature points. Although computation time varied by almost a factor of 
three, there was little effect on the accuracy, indicating that, for this trajec- 
tory at least, 12 points were sufficient. Thus the computation times shown 
in Table 9 might be reduced with no significant loss of accuracy by reducing 
the number, of. quadrature points. 


Finally, the programs used to evaluate the asymptotic solutions are still 
in the development stage and no real attempt has been made to optimize the 
computation time. The. programs make many calculations which do not 

I " , V 

affect .the ..solution but which we re useful during checkout. It is estimated 
that by eliminating these unnecessary calculations and making other changes 
within. the program, the computation times could be reduced from 25 to 50 
percent!. ...... ... ., , .. • ... - • 

> Vi' • ’ i A* . • r ' , * A , 

4. 6 DISCUSSION OF NUMERICAL RESULTS 

The numerical results presented he, re consist of a limited number of 
example s. from which various trends in accuracy can be observed. These 
trends are.fe.lt to be. representative of. the accuracies which result from the 
different boundary value, solutions, but it may be that a more comprehensive 
numerical analysis might reveal. further information. 
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The examples for which data is presented include those which were used as 
test cases in checkout" of the computer programs for each solution. This 
process went through the following sequence of solutions: 

1. Linear first-order midcourse interplanetary 

2. Linear first-order interplanetary 

3 . ' Linear second-order midcourse interplanetary 

4. Linear second-order interplanetary 

-5. Modified linear second-order interplanetary 

6. Linear second-order Earth-to-moon 

7. Nonlinear second-order Earth-to-moon 

8. Linear second-order two-impulse moon-to-Earth. 

Since the analytical derivations of each solution preceded the coding and 
checkout by several months, the derivation of the final solution (8) was 
completed while the coding and checkout of 4 was taking place. Also; by the 
time results were being obtained from 6 and 7, the coding of 8 had been 
completed. Since many of the factors influencing the asymptotic solution 
were determined from the numerical analysis of 6 and 7 and since this 
analysis occurred late in the study period, there was not sufficient time to 
develop a more comprehensive numerical survey and make changes in the 
other solutions. 

The sequence of solutions did reveal two aspects of the asymptotic solution 
which need to be re-emphasized. The first is the need for a better nominal 
or zeroth-order outer solution in order to decrease the deviation between 
it and the first order solution. An attempt was made to improve the nominal 
in the linear midcourse interplanetary solution. The new nominal did ndt go 
through the center of the target body and thus (2. 2-4) was violated. This 
caused singularities in the definite integrals and therefore an invalid solution 
This was corrected by introducing a fictitious body with the mass of the ' 
target body but with a slightly different position in order to make the new 
nominal satisfy (2. 2-4). This hew solution did give slightly better 
results for Cases 301 and 302 and eventually led to the inclusion of the 
nonlinear versions of each solution. [The nonlinear versions were first 
discussed by Carlson (Reference 2) but were not emphasized strongly.] 
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The second characteristic which was revealed was that the form of the 
fundamental solution can affect the accuracy. It was mentioned in Sub- 
section 2.6 that the form chosen was not unique and was slightly different 
than the results of Carlson (Reference 2). A more comprehensive numerical 
survey should include different versions based on these differences in the 
fundamental solution. 

A general comment regarding all of the solutions is that attempts .to improve 
accuracy, by whatever means, sometimes results in an improvement in one 
boundary condition and a degradation in another. This makes it somewhat 
difficult to assess the value of the supposed improvement without studying a 
large number of cases. 

It would appear from the results obtained in the preceding sections and from 
the comments of this section that a more comprehensive numerical survey 
of all of the solutions is needed. Such a survey should be based essentially 
on the solutions as they have been derived but including all versions of each 
solution and those changes, such as the construction of a composite-type 
solution, which have been specifically mentioned previously. 
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Section 5 

CONCLUSIONS AND RECOMMENDATIONS 
5. 1 CONCLUSIONS 

The analytical solutions presented in this report are a result of combining 
and extending the best ideas put forth in previous works on this subject. 

The solution of Reference 1 stressed the mathematical rigor and explicitness 
of results which evolved from the earlier work of Lagerstron and Kevorkian. 

f 

Reference 2, while also stressing a certain degree of explicitness, was based 
on a more general solution first derived by Breakwell and Perko. The 
concepts of mathematical rigor, explicitness, and generality have been 
combined in this study to generate a number of asymptotic solutions for both 
lunar and interplanetary applications. The basic solutions are the second 
order outer and inner solutions presented in Subsections 2. 1 and 2. 2. These 
solutions are themselves useful for studying motion along perturbed ellipses 
and hyperbolas, respectively; however, their main purpose in this study lies 
in the formulation of general solutions for trajectories which have both 
elliptic and hyperbolic behavior, with a continuous transition from one to the 
other. 

The fundamental solution presented in Subsection 2. 5 represents the result of 
matching the second-order outer and inner solutions in the overlap domain 
(i. e, , the region of transition from elliptic to hyperbolic motion). This 
solution is valid for any number of perturbing bodies and represents the first 
matched asymptotic solution which has been carried to second order. 

Several boundary value solutions have been formulated directly from the 
fundamental solution. These include the Earth-to-moon solutions (Sub- 
sections 3.2 and 3.3) and the interplanetary midcourse solution (Subsec- 
tion 3. 6). These solutions are second-order approximations to the exact 
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solution and satisfy the prescribed boundary conditions (both initial and 
terminal) exactly without need for any iterative techniques (other than in the 
Lambert solutions). 

Two fundamental solutions have been combined to formulate the interplanetary 
solution of Subsection 3. 7. Like the other boundary value solutions it is 
N-body, second-order, and non-iterative. It is valid for all trajectories 
originating close to one planet and terminating close to another. 

k. * • i I , 

A modified version of the fundamental solution has been used to formulate 
the one-impulse moon-to-Earth solution (Subsection 3. 4). In addition, 
separate versions of the second-order inner and outer solutions’ have been 
used directly (no fundamental solution) to formulate the two-impulse moon- 
to-Earth solution (Subsection 3. 5). ' 

Although applicable to a number of boundary value problems, all of the 
solutions presented are similar in that they are represented analytically 
in the form of asymptotic expansions. Also, the solutions have been derived 
so that most of the mathematical operations required for numerical evalua- 
tion are common to each formulation and all of the solutions can be program- 
med as subroutines of a master trajectory program. 

The primary purpose of this study was the development of analytical 
boundary value solutions. Nevertheless, Section 4 contains sufficient 
comparisons with numerical integration to demonstrate the accuracy and 
speed of the asymptotic solution in typical applications. For each solution 
evaluated, the results showed a significant improvement in accuracy over 
standard conic approximations. In one application (i. e. , interplanetary 
midcourse the accuracy is comparable to that of numerical integration. 

(In fact, the results of the asymptotic solution pointed out a deficiency in the 
numerical integration program with which it was being compared. Introduc- 
tion of more stringent internal accuracy requirements in the numerical 1 
integration program rerhoved the deficiency arid brought the numerical 
results into agreement with those of the asymptotic solution, ) . 
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The computation times for all of the solutions are significantly less than 
those for numerical integration and in some instances are less than the 
computation times of the fastest numerical approximation techniques. This 
computational speed makes the asymptotic solution very attractive for 
boundary value problems since all the numerical techniques, exact or 
approximate, require iterative methods to converge to a solution. 

The best overall application for the asymptotic solution appears to be the 
midcourse application, particularly for interplanetary trajectories where 
p is small. In this application the solution is convergent to second order, 
i. e. , the second-order error is less than first order. The solution is also 
less subject to sensitivity problems which arise with lunar trajectories 
originating close to perigee or close to the moon and with interplanetary 
trajectories originating close to the launch planet. The lack of uniform 
convergence of the second order solution and the problem of sensitivity are 
closely associated and both have a minimum effect in the midcourse 
solutions. 

o 

In other applications certain boundary conditions may cause the second- 
order error to be larger than first order. From an examination of all 
the data generated (more than has been presented in Section 4) it appears 
that the cause of the second-order divergence, when it occurs, is not a 
deficiency in the. second-order terms themselves, but rather the fact that 
the zeroth-order solution is a poor approximation to thej actual trajectory. 

This problem is discussed several times in Section 4 and again in the next 
section. , 

I 

The general conclusions are: (1) the second-order asymptotic solution does 
converge (with certain exceptions), (2) it is significantly more accurate than 
conic approximations., and (3) it is much faster than numerical integration. 

5. 2 RECOMMENDATIONS FOR FURTHER STUDY , 

In order to formulate a second-order asymptotic solution with more uniform 
accuracy, ft is recommended that any further studies include the following: 

A. A more comprehensive numerical survey which would include 

complete families of trajectories generated by varying each boundary 
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condition through a wide range of values and, in the interplanetary 
cases, trajectories to several different planets. Such a survey 
would better define those conditions under which the different 
second-order solutions diverge and would indicate exactly where a 
better nominal or zeroth-order solution is required. 

B. An analysis to determine what modifications can be made to the 

zeroth-order outer solution to make it a better approximation to the 
actual trajectory in those cases where divergence of the second-order 
solution does occur. This problem is discussed in some detail in 
, Section 4, but it should be re-emphasized that (2.2-4), which is an 
essential condition to successful matching of the outer and inner 
solutions, restricts the zeroth-order solution in such a way as to 
have an adverse effect on the second-order accuracy. It may be 
possible, however, to construct a piecewise continuous zeroth^ 
order solution from two-body arcs passing through a sequence of 
points obtained from the first-order solution. Only the arc closest 
to the target (or launch) body would have to satisfy (2.2-4) yet the 
total sequence, since it is obtained from the first-order solution, 
would be a much better approximation than the original nominal 
solution. Surh an approach would cause some increase in computa- 
tion time but should result in a uniformly valid solution for all 
applications, that is, a solution free of second-order divergence. 

Although A and B are related, they could be carried out independently. If 
such is the case, then B should be given the higher priority since it would 
produce the more useful results. 
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